(navigation image)
Home American Libraries | Canadian Libraries | Universal Library | Community Texts | Project Gutenberg | Children's Library | Biodiversity Heritage Library | Additional Collections
Search: Advanced Search
Anonymous User (login or join us)
Upload
See other formats

Full text of "Multiscale spline derivative-based transform for image fusion and enhancement"

A MILTISCALK SPLINE DERIVATIVE-BASF!) TRANSFORM 
FOR IMAGE FUSION AND ENHANCEMENT 



By 
IZTOK KOREN 



A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL 

OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT 

OF THE REQUIREMENTS FOR THE DEGREE OF 

DOCTOR OF PHILOSOPHY 



UNIVERSITY OF FLORIDA 



1996 



Copyright 1996 



by 



Iztok Koren 



Mojim starsem 



ACKNOWLEDGEMENTS 



I would like to thank my mentors Dr. Fred J. Taylor and Dr. Andrew F. Laine 
for their help and support. They gave me the freedom to explore research directions 
that interested me most and provided me with the excellent work environment. 

I am grateful to Dr. Jose C. Principe, Dr. John M. M. Anderson, and Dr. 
Kermit Sigmon for serving on my committee and for their comments on my disser- 
tation. 

I would also like to thank all the members of the High Speed Digital Architec- 
ture Laboratory and the Image Processing Research Group for their friendship and 
help. 

Last but not least, I want to thank my family for their support and under- 
standing, and for bearing with me when I was unbearable. 



IV 



TABLE OF CONTENTS 



ACKNOWLEDGEMENTS iv 

LIST OF TABLES vii 

LIST OF FIGURES viii 

ABSTRACT ix 

CHAPTERS 

1 INTRODUCTION 1 

1 . 1 Shortcomings of Traditional Methods of Wavelet Analysis .... 1 

1.2 Thesis Overview 4 

1.3 Notation 5 

2 SIGNAL PROCESSING USING CENTRAL B-SPLINES 8 

2.1 Central B-Splines: Definition and Properties 8 

2.2 B-Spline Signal Interpolations 9 

2.3 B-Spline Signal Approximations 14 

3 TRANSFORM IN ONE DIMENSION 17 

3.1 1-D Discrete Dyadic Wavelet Transform Revisited 17 

3.2 Remarks 28 

4 TRANSFORMS IN TWO DIMENSIONS 30 

4.1 2-D Discrete Dyadic Wavelet Transform Revisited 30 

4.2 Steerable Functions 36 

4.3 Steerable Dyadic Wavelet Transform 38 

4.4 Multiscale Spline Derivative-Based Transform 40 

5 IMPLEMENTATION ISSUES 49 

5.1 Finite Impulse Response Filters 49 

5.2 Infinite Impulse Response Filters 55 

6 APPLICATIONS 59 

6.1 Image Fusion 59 



6.2 Image Enhancement 61 

7 CONCLUSION 66 

REFERENCES 68 

BIOGRAPHICAL SKETCH 72 



VI 



LIST OF TABLES 



2.1 Transfer functions of direct B-spline filters for orders from to 9. . . 12 

3.1 Impulse responses h(n), g(n), l(n), and k(n) for p = and d € {1,2,3}. 25 

3.2 Impulse responses h(n), g(n), /(??), and k(n) for p=l and d € {1,2.3}. 25 

3.3 Impulse responses h(n), g(n), /(n), and k(n) for p — 2 and c/ € {1,2,3}. 26 

4.1 Impulse responses tx(n) for p € {0,1,2} 32 



VII 



LIST OF FIGURES 



1.1 Discrete wavelet transform of two signals translated to each other. . . 2 

1.2 Discrete wavelet transform and "algorithme a trous." 3 

2.1 Spline functions p (x) and r) p (x) for p € {0,1,2,3,4} 10 

2.2 Fourier transforms 3 p (ui) and rj p (u>) for p € {0,1,2,3,4} 13 

2.3 Splines /5 p (x) and /3 p (w) for p € {0,1,2,3,4} 16 

3.1 Filter bank implementation of a 1-D discrete dyadic wavelet transform. 21 

3.2 Comparison of two discrete implementations using sinc(.r) as an input. 23 

3.3 Wavelets r(.r) = d \j (x) for p € {0. 1,2} and d G {1.2.3} 27 

4.1 Filter bank implementation of a 2-D discrete dyadic wavelet transform. 35 

4.2 Filter bank implementation of a spline derivative-based transform. . . 48 

6.1 Comparison of image fusion results 62 

6.2 Comparison of image enhancement results 64 

6.3 The magnitude of filter coefficients for G$(uJ x ,ujy) 65 



vin 



Abstract of Dissertation Presented to the Graduate School 

of the University of Florida in Partial Fulfillment of the 

Requirements for the Degree of Doctor of Philosophy 



A MULTISCALE SPLINE DERIVATIVE-BASED TRANSFORM 
FOR IMAGE FUSION AND ENHANCEMENT 

By 

Iztok Koren 

December 1996 

Chairman: Dr. Fred J. Taylor 

Cochairman: Dr. Andrew F. Laine 

Major Department: Electrical and Computer Engineering 

Analyzing images along distinct scales is advantageous in many image pro- 
cessing and computational vision tasks. Traditional wavelet transforms have become 
a popular tool for multiscale image analysis and image compression suffer from arti- 
facts, lack shift and rotation invariance, and may introduce aliasing and anisotropies. 
We propose a highly redundant representation that eliminates these artifacts. 

A new transform is constructed as an extension of the discrete dyadic wavelet 
transform with a wavelet being the first derivative of a central B-spline. We begin 
the construction in one dimension by providing an efficient initialization procedure 
for the computation of the discrete dyadic wavelet transform and by extending the 
transform to encompass higher order derivatives. 

This modified one-dimensional discrete dyadic wavelet transform will serve 
as a basis for deriving a steerable two-dimensional transform. Such a transform is 



IX 



advantageous in that it does not introduce any of the artifacts described above. How- 
ever, exact reconstruction is problematic via implementation in the spatial domain. 
We solve this problem by devising a spline-based approximation. The resulting algo- 
rithm for computing a multiscale spline derivative-based transform is implemented 
efficiently as a filter bank consisting of separable filters alone. For both finite and 
infinite impulse response filters present in the filter bank, we specify implementa- 
tion details of a realization that alleviates the boundary effects caused by the use of 
circular convolution. 

Finally, we use the derived transform for multiscale image fusion and enhance- 
ment. We show empirically that the transform does not introduce artifacts commonly 
reported for traditional wavelet-based implementations and provides more flexibility 
for different fusion and enhancement strategies by enabling orientation analysis. 



CHAPTER 1 
INTRODUCTION 



Analyzing images across multiple scales and resolution has become a power- 
ful tool for solving compelling problems in computational vision, image processing, 
and pattern recognition. Wavelet theory encompasses multiscale and multiresolution 
representations, such as subband filtering [40], image pyramids [4], and scale space 
filtering [48], into a unified mathematical framework. In the area of image processing, 
there remain few research areas to which wavelet analysis has not been applied. For 
example, problems in image compression, denoising, restoration, enhancement, reg- 
istration, fusion, segmentation, and analysis, have all been approached with distinct 
kinds of wavelet processing. 

Though ubiquitous, wavelet analysis is not without problems of its own. Lack 
of translation invariance, one of the major problems of the wavelet transform [10], is 
in multiple dimensions accompanied with lack of rotation invariance. 

1.1 Shortcomings of Traditional Methods of Wavelet Analysis 

A wavelet transform in its most commonly used orthogonal or biorthogonal 
forms is not translation and rotation-invariant. By translation-invariant transform, 
we mean a transform that commutes with a translation operator. Since we will 
deal primarily with discrete transforms in this work, we constrain the translation 
parameter to integer multiples of a sampling period. 

Lack of translation invariance of the discrete wavelet transform is illustrated 
in Figure 1.1. Here, we can clearly see how a translation of the input signal by one 





■ 




4 ill* 



II 

, ! 


i 4 in* 

6 





9 I r 




(a) 



(b) 



Figure 1.1: (a) Original signal and (b) signal translated one sample to the left with its 
discrete wavelet transform coefficients shown across dyadic scales 2 m , m G {1,2,3,4}. 



sample results in a completely different set of transform coefficients (orthogonal 
wavelet DAUB4 1 [10] was used in this experiment). 

Noninvariance under translations of an orthogonal and biorthogonal wavelet 
transform is due to lower sampling density at coarser scales. 2 A straightforward way 
of dealing with this problem is to construct a redundant transform by using the same 
sampling frequency for the input signal and all scales of the transform. A filter bank 
implementation of such a transform, called "algorithme a trous" [15], is based upon 
the fact that downsampling followed by filtering is equivalent to filtering with the 
upsampled filter before the downsampling, as shown in Figure 1.2. 



G(z) 2t 



H(z) 2t 



G(z) 2t 



G(z) 



G(z) 2t 



H(z) 



G(r) 



H(z) 2t 



(a) 



G(z 4 ) 



H(r) 



H(z') 



(b) 



Figure 1.2: Filter bank implementation for (a) a discrete wavelet transform and (b) 
"algorithme a trous" decompositions for three levels of analysis. 



Lack of rotation invariance is another shortcoming of traditional (i.e.. orthog- 
onal and biorthogonal) wavelet techniques. In defining rotation invariance, we are a 
bit less strict than with translation invariance. We do not require that the transform 
commutes with a rotation operator here. Even in the case of a simple filtering, this 
would limit us to circularly symmetric filters only. Our requirement for analysis is a 



lr rhe number in DAUB4 refers to twice the order of the wavelet (i.e.. two in this case). 
In practice, since analysis is performed over a finite range of scales, a discrete wavelet transform 
is translation-invariant by translations determined by the coarsest scale (e.g., sixteen samples for 
the analysis from Figure 1.1) [10]. 



transform that enables rotation-invariant processing. As an example of such a trans- 
form, let us consider filtering with the first derivative of a two-dimensional Gaussian 
probability density function in two directions, specifically, along x and along y-axis. 
By linearly combining the results of filtering in these two directions, filtering with the 
first derivative of a Gaussian in any direction can be computed. This fact was used 
by Canny [6] for edge detection. A determined edge direction rotates as an input 
image is rotated. 

After choosing the fundamental properties of the transform, one must decide 
upon the basis functions to be applied. For our studies, we selected basis functions 
that well approximated derivatives of a Gaussian, because (1) the Gaussian proba- 
bility density function is optimally concentrated in both time and frequency domain, 
and thus suitable for time-frequency analysis, (2) higher order derivatives of a Gaus- 
sian can be, similar to the first derivative, used for rotation-invariant processing [12], 
and (3) the Gaussian function generates a causal (in a sense that a coarse scale de- 
pends exclusively on the previous finer scale) scale space [2]. The last property makes 
possible scale-space "tracking" of emergent features. 

1.2 Thesis Overview 

We wish to construct an invertible, translation and rotation-invariant two- 
dimensional transform that decomposes an image using a set of basis functions which 
are approximations to derivatives of the Gaussian probability density function. 

To obtain a translation-invariant representation, we will use "algorithme a 
trous," or more precisely, we will extend the discrete dyadic wavelet transform [28], 
which uses "algorithme a trous." Derivatives of a Gaussian will be approximated by 
wavelets which are derivatives of central B-spline functions. In Chapter 2, we review 
the basics of B-spline processing necessary to derive transforms in later chapters. 

Next, Chapter 3 examines a dyadic wavelet transform in one dimension. We 
present extensions of the discrete dyadic wavelet transform to higher order derivatives 



and to even order spline functions. To compute a discrete transform from a contin- 
uous one, the discrete computation must be properly initialized. We devise a new 
initialization procedure that is shown to be more accurate than the one suggested by 
Mallat and Zhong [28]. After the derivation of the transform, we point out relevant 
connections to scale-space filtering and reconstruction from edges. 

In case of the first derivative of a Gaussian, a rot at ion- in variant transform 
can be obtained by a tensor product extension of the one-dimensional transform 
to two dimensions. For second order derivatives, a tensor product approach can 
approximate Laplacian of a Gaussian, but cannot perform orientation analysis (due 
to its isotropic nature). In Chapter 4, we first examine the above two cases, and 
then develop a new multiscale spline derivative-based transform, which in addition 
enables rotation-invariant processing for derivatives of orders higher than two. 

Given that the developed transform is highly redundant, an efficient imple- 
mentation is extremely advantageous. This issue is addressed in Chapter 5. When 
filtering of a signal is performed by circular convolution, boundary effects may result. 
A standard way of dealing with this problem is by mirror-extending the input signal 
to the filter. We combine such an extension with symmetry/antisymmetry of the 
filters, to achieve savings in both computation time and memory. 

Finally, in Chapter 6, we present two sample applications: image fusion and 
enhancement. We demonstrate that the new transform does not cause artifacts com- 
monly associated with traditional wavelet techniques and is an attractive analytic 
tool for designing novel fusion and enhancement strategies. 

1.3 Notation 

We use symbols TV, Z, and R for the sets of naturals, integers, and reals, 
respectively. 

L 2 {R) and L 2 (R ) denote the Hilbert spaces of measurable, square-integrable 
functions f(x) and f(x,y), respectively. 



The inner product of two functions f(x) £ L 2 (R) and g(x) £ L 2 (R) is given 

by 

(f(x),g(x)) = f°° f(x)g(x)dx. 

J — oo 

The norm of a function /(.r) £ L 2 (R) is defined as 



11/11 = JJ~Jf(x)\*dx. 

The convolution of functions /(x) £ L 2 {R) and </(x) £ L 2 (R) is computed as 
/*</(-*•)= /" f(t)g(x-t)dt 

J — 'OO 

and the convolution of two functions f(x,y) £ L 2 {R 2 ) and g{x,y) £ L 2 (R 2 ) equals 

f*g(x,y)= / f(tx,ty)g(x-t x ,y-t y )dt x dt y . 

J—oo J—oo 

The Fourier transform of a function /(.r) £ L 2 {R) is defined as 

/(w)= / f(x)e-'»*dx, 

J — oo 

and the Fourier transform of a function f(x,y) £ L 2 {R 2 ) is equal to 
/(u/„w y ) = r r f(x,y)e-^ + ^dxdy. 

J —'OO J — oo 

l 2 (Z) and 1 2 {Z 2 ) stand for the spaces of square-summable discrete signals f(n) 
and f(n r ,riy). respectively. 

The --transform of a discrete signal f(n) £ 1 2 {Z) is defined as 

F(z)= £ f(n)z-\ 

n=—'Oo 

The convolution of discrete signals f(n) £ P(Z) and g(n) £ / 2 (Z) is equal to 

oo 

f*g(n)= ^2 f( m )9(n-m), 

m=— oc 

and the convolution of discrete signals f(n x ,n y ) £ l 2 (Z 2 ) and g{n x ,n y ) £ 1 2 {Z 2 ) is 
given by 



oo oo 



/ * g(n x ,n y )= Yl Yl f( m *> m y)g{ n r- m x , n y - m y ) . 



m x =—oo m-u=—oo 



The Fourier transform of a discrete signal f(n) £ l 2 (Z) is equal to the 
transform evaluated on the unit circle 

CO 

F{u>)= £ /(»>-^ n , 

n= — oo 

and the Fourier transform of a discrete signal f(n x ,n y ) £ l 2 (Z ) is defined as 
FKw s )= £ £ /(n„n lr ) C -^- B - +< -» B »>. 
For later use, we define the following functions: 

1. the unit impulse function 

r . . 1 for x = 

I I) otherwise, 

2. the unit step function 

( 1 for x > 
U{X) := for x < 0. 



3. the rectangular function 



4. the sine function 



rect(x) := < n , ? 

for x > i 



• / v sin(7r.r) 
smc(.r) := , and 

1TX 



5. the unit impulse sequence 

£( \ J 1 ^° r ?? = 

I otherwise, 
where x £ jR and n £ Z . 



CHAPTER 2 
SIGNAL PROCESSING USING CENTRAL B-SPLINES 



In this chapter we briefly review fundamentals of spline processing needed for 
derivations in Chapters 3 and 4. 

2.1 Central B-Splines: Definition and Properties 

Given real numbers — oo < x < x\ < ,r 2 < . . . < x m < x m+ i < oo, a function 
on the interval [.ro,.r m +i] is called a spline function of order p with the knot (i.e., 
grid point) sequence Xi, x-2. ... .r m . if it is (1) a polynomial of degree p or less in each 
interval [a?,-, x i+ i]. i = 0, 1, . . . m, and (2) continuous in its derivatives up to the order 
p-1 on the interval [.r ,.r m +i] (i.e., C p ~ 1 [x ,Xm+i})- 

Here, we will concentrate primarily on basis splines (B-splines), or more pre- 
cisely, central B-splines having knots at i £ Z for p odd and at i + | for p even [34]. 
Central B-splines of order p (with p+1 knots) are defined as 

Figure 2.1 shows 3 p (x) and their Fourier transforms (5 p {u)) for p £ {0, 1, 2, 3, 4}. 

A family of functions {/3 p (x — ???)} % forms a basis of S p , a space of order 
p spline functions with knots at i for p odd and at / + | for p even (?' G Z) [34. 35]. 
Except for p = 0, the basis functions {/? p (x — m)} are not orthogonal. 

Let us list some properties of functions j3 p [x) [3, 34]: 

1. fip{x) are nonnegative functions with a support of length p+1, 



p-\-l times 

/ * V 

2. /# p (.r) = fio * /?o * • • • * A)(#) 5 where "*" denotes the convolution operator, or, 
equivalently, in the Fourier domain: 

3. P P {x) = I ((*? + x) & P - X {x + \) + (*!-*) ^(x - §)), 

*-¥ = M j+ !)-M*-!)- 

Another interesting property of B-splines is the fact that they converge to a 
Gaussian as their order tends to infinity. Unser et al. [44] derived the Gaussian 
approximation 

i -4 

where a v — J 2 ^-- Ratio n7 ,n , for example, is already well below 0.2%. 

v i^ lift* 1 1 

Denoting by S\ a spline function space spanned by {fi p (x — m)} for p odd 
and by {/3 p (x — \ — rn)} for p even (subscript in Sj refers to the fact that spline 
functions have knots at integers), spline spaces form a nested sequence . . . C St, C 
. . . C S j C Sj C S^-i C . . . C S^-i C .... By orthogonalizing this basis functions 
a multiresolution analysis of L 2 (R), from which the Battle- Lemarie wavelet bases 
stem, can be built [10]. Here, we will not pursue constructions of orthogonal, semi- 
orthogonal, or biorthogonal spline-wavelets any further — reader looking for a detailed 
treatment of this subject may find a good starting point in books [8, 9]. 

2.2 B-Spline Signal Interpolations 

Unser et al. developed a fast digital filtering scheme for B-spline signal pro- 
cessing [43]. They defined a discrete B-spline of order p and expansion factor (spacing 
between knots) m as 

b p ,m(n) := ,8 P l—j , n,m € Z. 



10 



p-n 



-3 -2 -1 








I a | 



(b) 



Figure 2.1: Spline functions (a) j3 p (x) and (b) t} p (x) for p <E {0, 1,2,3,4}. 



II 

Henceforth, if the distance between knots equals one, we will write b p (n) instead of 

Interpolation of a discrete signal s(n) G l 2 (Z) by s p (x) G S p using central 
B-splines 



s(n) = s p {x) 



J2 c(i)/3 p (x-i) 



5 

x=n 



(2.1) 



can now be written as a convolution sum 

s(n) = c*b p (n). (2.2) 

If s(n) are samples of a function s(x) bandlimited to [— k. k] (i.e., the support 
of its Fourier transform S(uj) is in [— 7r,7r]), it can be shown that s p {x) — ► s(x) as 
p — > oo [1, 36]. 

In [13]. they refer to a linear operator by which B-spline coefficients c(n) can 
be obtained from samples s(n) as a "direct B-spline transform." Equation (2.2) 
therefore represents "indirect B-spline transform'' of a sequence {c(n)}. 

After taking the --transform of (2.2), the direct B-spline filters are found to be 
^p 1 ! 77 ) = Z~ 1 {[Bp( z )]~ 1 }- Since B-spline functions j3 p (x) are compactly supported, 
indirect B-spline filters b p (n) are finite impulse response (FIR) filters, while direct 
B-spline filters ^(n) are infinite impulse response (IIR) filters. Aldroubi et al. [1] 
showed that IIR filters b~ l {n) are stable (i.e., the region of convergence of B~ l {z) 
includes the unit circle [30] ) for any order p. Note that both indirect and direct B- 
spline filters are symmetric, which follows from the fact that central B-splines fi p {x) 
are symmetric. 

Table 2.1 shows the --transforms of direct B-spline filters for the first ten 
orders. We postpone the discussion on implementation details of B-spline filters until 
Chapter 5. 



L2 



Table 2.1: Transfer functions of direct B-spline filters for orders from to 9. 



vo 



8 



6 
z+4+z -1 

384 



~? +76J+230+76J- 1 +z~ 2 

120 

j 2 +26-+66+26;- 1 +.— 2 

46080 



'4-722; 2 + 10543j4-23548+10.543j- 1 +T22j- 2 +j- 3 

5040 

j 3 + 120; 2 + 1191 i+2416+l 191 z~* +120z~ 2 +z- 3 

10321920 



: 4 +6552j 3 +331612- 2 +24852883+4675014+2485288;- 1 +331612z- 2 +6552^- 3 + ;- 4 



3(12880 



: 4 +502- 3 + 14608s 2 +88234z+156190+88234z- 1 +14608;- 2 +502:- 3 + -- 4 



Instead of using B-spline interpolation as given by (2.1) it is sometimes con- 
venient to express the interpolating function s p (x) in terms of discrete samples s(n) 



s p( x ) = Y, s ( l )Vpi x -i), 



(2.3) 



where r/ p (.r ) is the cardinal spline of order p. In the frequency domain, cardinal splines 
converge to an ideal lowpass filter with cutoff frequency tv (i.e., r] p (x) — ► sinc(.r)) as p 
tends to infinity [1. 45], which establishes the asymptotic equivalence with Shannon's 
sampling theorem [37]. 

Using (2.1) and (2.2) with (2.3) cardinal splines can be related to B-splines: 



V P (x)= E b;\i)P P (x-i). (2.4) 

i=—oo 

Cardinal splines T) p (x) and rj p {w) for p G {0, 1,2,3,4} are shown in Figure 2.2. 



13 




-2s 2s 




-2» 2s 




-8s -6s ^ -2s 2s 4* 6a 8n 




-6s -4s -2s 2s 4™ 6* 8s 




(a) 




-6s -4* -21 2s 4: 6x 8* 




-8s -6s -4s -2 7 2* 4s 6s 6* 




-8s -6s -4* -2* 2s 4s 5j 




-6s -6s -4s -2s 2s 4s 6s 8s 




-8s -6s -4s -2s 2. 4s 6s 8s 



(b) 



Figure 2.2: Fourier transforms of spline functions (a) 3 p {lj) and (b) rjJu) for p £ 

{0,1,2,3,4}. 



II 

2.3 B-Spline Signal Approximations 

Central B-splines are also simple to use when the goal is function approxima- 
tion. Least-squares B-spline approximation of s(x) £ L 2 {R) is achieved by computing 
the orthogonal projection of this function onto S p . We have 

~s(x) = s p (x) = f] d(i)/3 p (x-i) (2.5) 

i=—oo 

with 

d(i) = (s(x)J p {x-i)), (2.6) 

where 

°° 

&(*) = £ a&iCO&C* - 

!= — OO 

O 

is the dual spline of order p [45]. Spline functions /3 p (x) and /? p (:r) satisfy the biorthog- 
onality condition 

o 

((3 p (x - m), i3 p (x - n)) = 6(m - n), m, n e Z. 



(Note that, since both /3 p (x) and ft p (x) form a basis of 5 P , they can be interchanged 

in (2.5) and (2.6).) 

O o 

Figure 2.3 shows functions f3 p (x) and their Fourier transforms 3 p (uj) for p 6 
{0,1,2,3,4}. 

An interesting alternative to the minimum Z 2 -norm (i.e., least-squares) ap- 
proximation of a signal is obtained by computing oblique instead of orthogonal pro- 
jection of the signal onto the spline function space. Unser and Aldroubi proposed an 
independent specification of the sampling and approximation spaces [42]: a linear op- 
erator maps coefficients of the input signal expansion over sampling space basis into 
the coefficients of the approximation space basis expansion from which the signal's 
projection onto the approximation space is recovered. Constraining the entire system 
to be linear, shift-invariant for integer translations, and consistent (i.e., the system 
acts as an identity operator for functions that belong to the approximation space), 



l.-J 

the obtained solution for the signal approximation is the projection of a signal onto 
the approximation space perpendicular to the sampling space [42]. Analogously to 
(2.5) and (2.6) this projection can be expressed as 

sr(x)= f] a{i)0 r (x-i) (2.7) 

i=— oo 

with 

a(i) = q£ *(s(x),p 3 (x-i)), 

where q^} is the convolution inverse of the cross-correlation sequence 

?«(») = {P$(x - i),Pr(x)) = b s+r+1 {i). 

When the sampling space S s and the approximation space S r are identical 
(i.e., 5 = r), an orthogonal projection given by (2.5) and (2.6) results. (Note that 
the described oblique projection is not restricted to spline function spaces — the only 
requirement is that both sampling space basis and approximation space basis are 
Riesz bases of the corresponding function spaces [42].) 

Signal approximation (2.7) is particularly attractive in situations where the 
sampling space is given a priori (e.g., by the impulse response of the acquisition device 
[42]) or when such an approximation is close to the optimal least-squares solution but 
simpler to implement than the orthogonal projection (e.g., [47]). 



I(i 




-15 -to -5 



-15 -tO -5 



S 10 IS 




5 10 




5 tO 15 




10 15 





0-1 

s 

■ 


-8* -6« -4* -2i ir 


P-2 

. _J 


, 





-At 2* 4» &> Si 



-<*I -2* 2i 4i 6i b 



5 10 15 




(a) 

o 

Figure 2.3: Spline functions (a) /3 p (x) and (b) their respective Fourier transforms 
4(«)forp€ {0,1,2,3,4}. 



CHAPTER 3 
TRANSFORM IN ONE DIMENSION 



In this chapter, we will augment the one-dimensional discrete dyadic wavelet 
transform [28] to obtain a strong foundation for derivations of two-dimensional trans- 
forms in Chapter 4. 

3.1 1-D Discrete Dyadic Wavelet Transform Revisited 

A discrete wavelet transform is obtained from a continuous representation 
by discretizing dilation and translation parameters such that the resulting set of 
wavelets constitutes a frame. The dilation parameter is typically discretized by an 
exponential sampling with a fixed dilation step and the translation parameter by 
integer multiples of a fixed step [10]. Unfortunately, the resulting transform is variant 
under translations, a property which makes it less attractive for image analysis (cf. 
Section 1.1). 

As we have already mentioned in Chapter 1, sampling the translation param- 
eter with the same sampling period as the input function to the transform results in 
a translation-invariant, but redundant representation. The dyadic wavelet transform 
proposed by Mallat and Zhong [28] is one such representation. Let us begin with a 
brief review of properties of the dyadic wavelet transform as described in [28], but 
included here for completeness. 

The dyadic wavelet transform of a function s(x) € L 2 (R) is defined as a 
sequence of functions 

nM^)} rae z- (3-1) 



IT 



IN 

where W m s(x) = s * i^ m (x), and m (.T) = 2 -m ^(2 -m x) is a wavelet tp(x) expanded 
by a dilation parameter (or scale) 2 m . Note the use of convolution instead of an inner 
product. 

To ensure coverage of the frequency axis the requirement on the Fourier trans- 
form of V'm(-c) is the existence of A\ > and B\ < oo such that 

oo 

A, < £ |0(2 m cu)| 2 < B x 

m=— oo 

is satisfied almost everywhere. The constraint on the Fourier transform of the 
(nonunique) reconstructing function \(x) is 

m=— oo 

A function s(x) can then be completely reconstructed from its dyadic wavelet trans- 
form using the identity 

oo 

S ( X ) = Yl W ^ S * Xm(x), 

m=— oo 

where \ m {x) = 2- m \(2- m .r). 

In numerical applications, processing is performed on discrete rather than 
continuous functions. When the function to be transformed is in the discrete form, 
the scale 2 m can no longer vary over all m £ Z. Finite sampling rate prohibits the 
scale from being arbitrarily small, while computational resources restrict the use of 
an arbitrarily large scale. Let the finest scale be normalized to 1 and the coarsest 
scale set to 2 . 

The smoothing of a function s(x) G L 2 {R) is defined as 

S m s(x) = s*<f> m (x), 

where <f> m (x) = 2~ m (p(2~ m x) with m G Z, and <j>(x) is a smoothing function (i.e., its 
integral is equal to 1 and <f>(x) -+ as Ixl — > oo). 



1!) 

In [28], a real smoothing function (f>(x) was selected, whose Fourier transform 
satisfied 

|^)| 2 = f)^(2 m w)x(2 m o;). (3.2) 

m=l 

In addition, it was shown that any discrete function of finite energy s(n) £ 1 2 {Z) 
can be written as the uniform sampling of some function smoothed at scale 1, i.e., 
s(n) = Sof(n), where f(x) G L 2 (R) is not unique. Thus, the discrete dyadic wavelet 
transform of s(n) for any coarse scale 2 M was defined as a sequence of discrete func- 
tions 

{S M f(n + a), {W m f(n + s)} m£[hM] } neZ , 

where s is a ir{x) dependent sampling shift. 

The above initialization s(n) - S f(n) is rather standard in the discrete 
wavelet transform computation [10]. although it yields correct results (i.e., the dis- 
crete wavelet transform is equal to the samples of its continuous counterpart) only 
when s(n) = Sqs(ti). Here, we will concentrate on wavelets which are derivatives of 
spline functions and this will lead us to a simple initialization procedure [46] that 
alleviates the above problem. 

For a certain choice of wavelets the discrete dyadic wavelet transform can 
be implemented within a fast hierarchical digital filtering scheme. Next, we shall 
summarize the relations between filters, wavelets, and smoothing functions. 

First, let us introduce a real smoothing function <p(x) such that (3.2) can be 
rewritten as 1 

oo 

#(w)0M= £^(2"V) x (2 m a;), (3.3) 

m=0 

and let us set 6{x) = /3 p (x) (i.e., we restrict ourselves to wavelets which are spline 
functions). 



^ote that the sum index determines the range of scales of the discrete transform: using (3.2) 
we have ip(2u>) and \(2u) at the finest scale of the transform, while for (3.3) we get ip(u) and \(u). 



20 

Computing (3.3) for the finest two scales shows that 

0(w) x(u) = A(«)^(«) - &(2w) 0(2w). (3.4) 

/5 p (2w) can be related to /i p (u;) by expressing /i p (2u>) as (cf. Proposition 1 of [46]) 

1 / sin(^) \ / sin UJ 



p (2w) 



2P+1 \™(l). 



and using ELo^'"^ = !i ^lf 1 e4 ^ ): 

4,(2w)=(oob(|)) P+1 A(«). (3.5) 

(Note that a relation similar to (3.5) can be derived for integer scales provided that 
the dilation parameter and the order p are not both even [46].) 

Let F(m) be a digital filter frequency response and let F s (uj) = e JUJS F(io). 

If we choose 

0(w) = G_ # (w)4p(«), (3.6) 

tp(2u)=L,(u)<p(u) t (3.7) 

x(u) = K.(u)<p(w), (3.8) 



and 



H(u:) = e^(cos(f)Y'\ (3 " 9) 



where s G {0. |} is a filter dependent sampling shift needed for g(n). /(/?). k(n), and 
h(n) to be FIR filters, and insert Equations (3.5)-(3.9) into (3.4), we observe the 
relation between frequency responses of the filters 

G(u)K(u) + H(u)L(u) = 1. (3.10) 

Similar to orthogonal and biorthogonal discrete wavelet transforms, the dis- 
crete dyadic wavelet transform can be implemented within a hierarchical filtering 



21 



scheme. To derive such a digital filtering scheme, let us assume that s(iv) from (3.1 ) 
is bandlimited to [— 7r,7r]. Using Shannon's sampling theorem [37] and (3.6) in the 
definition of the dyadic wavelet transform (3.1) with m = shows 



II u 



/oo °° °° 
J2 s ( i )sinc( t - i ) ]T g_ s ( m ) 3 P ( x - t - m ) dt . 
-oo • _, _,_ 



By making use of the fact that the cardinal spline functions tend to the sine 
function as their order r approaches infinity, and employing (2.4) we can write 

Wo»H * S{u>)B?{u)P r (u>)$ p {a,)G-.{u,), 
or, by using (3.5) and (3.9), 

m — 1 

^{VU m .(.r)| x= J~.SV)^- 1 ^)i? P+r+ i^)G- 5 (2 m a;)n^- s (2 ! -'). (3.11) 

;=o 

Equation (3.11) entirely specifies the discrete dyadic wavelet transform de- 
composition, while the reconstruction follows from (3.4)-(3.9). Three levels of a filter 
bank implementation are shown in Figure 3.1. (Note that the initialization is the 
same as the one proposed in [46] and that except for the prefiltering and postfilter- 
ing, this scheme is implementing "algorithme a trous.") Noninteger shifts at scale 1 
for filters with s = | are rounded to the nearest integer. 



— ^o» BpWco) 



G>) 



- H«o) 



K,(a,) 



Q(2<0) 



H(2<o) 



G.(4<o) 



HJ4<:,) 



^(4(0) 



M4co) 



1^(20.) 



L s (2<o) 



M») 



^,«o); 8,(0.) 



Figure 3.1: Filter bank implementation of a one-dimensional discrete dyadic wavelet 
transform decomposition (left) and reconstruction (right) for three levels of analysis. 

We will now perform a simple experiment which will illustrate the difference 
between the implementation of the discrete dyadic wavelet transform as originally 



22 

proposed in [28] (i.e., without prefiltering and postfiltering) and the one from Figure 
3.1. 

Let s(x) = sinc(.r), p — 2, and g(n) = 28 (n + 1) — 28(n) (this particular choice 
for p and g(n) results in the same wavelet as was used by Mallat and collaborators 
[27, 28]). The dyadic wavelet transform of s(x) at a scale 2 m (3.1) in the frequency 
domain is then 

C*(w) = G. s {2 m u) fa(2 m u) rect (^) . (3.12) 

The Fourier transform of the discrete dyadic wavelet transform of s(n) = 8(n) at a 
scale 2 m using spline based initialization follows from (3.11) 

m-l 

F{W m s(n)} = B;\«:)B r+3 ^)G- s (2 m u,<) Y[H_ S (2'^). (3.13) 

i=0 

while the one using the algorithm from [28] equals 

m — 1 

F{W m s{n)} = G-,(2 m u>) J] H..(Tu). (3.14) 

i=0 

In Figure 3.2 a comparison of magnitudes of (3.13) and (3.14) versus (3.12) is 
shown: in Figure 3.2(a) magnitudes of (3.12) (solid) and (3.14) (dashed) are plotted 
for m € {0. 1,2.3}, while the dashed curves in 3.2(b) represent magnitudes of (3.13) 
with r = 5. 

By choosing the appropriate order r, (3.13) can approximate (3.12) in the 
interval [ — 7r, tt] arbitrarily good, while originally proposed (3.14) has troubles at finer 
scales. Mallat and Zhong [28] noticed that there was a problem with their discrete 
transform computation, and introduced a set of constants associated with the discrete 
transform coefficients at dyadic scales. They chose the values of constants such that 
the transform coefficient modulus maxima remained constant over all dyadic scales for 
a step edge input signal. In relation to Figure 3.2(a) this is equivalent to multiplying 
T{W m s(n)} by a distinct constant for each m. Clearly, this can improve over the 
situation depicted by Figure 3.2(a), but can still not compete with the described 
spline based initialization. 



23 




(a) 



(b) 



Figure 3.2: (a) Fourier transform magnitudes of the dyadic wavelet transform of 
s(x) = sinc(x) (solid) and the originally proposed discrete dyadic wavelet transform 
[28] of s(n) = S(n) (dashed), (b) Fourier transform magnitudes of the dyadic wavelet 
transform of s(x) (solid) and the discrete dyadic wavelet transform using quintic 
splines for interpolation of s(n) (dashed). 



24 

Next, we will choose filters in the filter bank implementation of the discrete 
dyadic wavelet transform. As already mentioned, we are interested in wavelets which 
are derivatives of spline functions. G(u>) in (3.6) is therefore the Fourier transform of 
the difference operator centered around —3 (cf. Property 4 in Section 2): 



G(w) = e*" 2; sin -J . (3.15) 



dmod2 



where d is the order of the derivative and the sampling shift for this filter is „ 

Since H(uj) was already given by (3.9), the remaining two filters to be deter- 
mined are L(u>) and K(ui). Both of them are (as is true for <p(x) and \{x)) nonunique. 
If we choose L(u>) such that we can express K[u) in terms of a finite geometric 
series having the smallest number of elements for an arbitrary p, we get 

L( U ) = e-^ y. (-i)- +l [ l ^ j \ ^ cos y j (3 16) 

and 

i d+1 < 

/ //.A \ dmod'2 / P / /, ,\ \ 2ra\ L 2 

' ,-Jujs ■ I w * 



where [x\ denotes the largest integer smaller than x, the sampling shift for L{u)) is 
the same as the one for H{u) (i.e., a= ^)^ od2 ) 1 and the sampling shift for K"(w) is 
the same as the one for G(u)). 

Note that Equations (3.9) and (3.15)-(3.17) work fine from the mathematical 
point of view, but in practice the reconstruction may become cumbersome when both 
p and d are large (the lengths of impulse responses ft(n), g(n), l(n), and k(n) are p+2, 
d+1, (p+l)(d-(d+l)mod2)+l, and pd+(p+l)(dmod2)+l, respectively, while for the 
frequency responses of the decomposition filters we observe that linip^oo \H- S {uj)\ = 
8 u (u + 2nir) and Yim d ->^(2j)- d |GL,(w)| = S u (u + {2n+l)ir) with n € Z). 

It is also worth noting that if(w) is a lowpass filter when p is even (i.e., the 
reconstruction function \{x) is a wavelet only for p odd). 



25 



Tables 3.1, 3.2, and 3.3 list impulse responses of the four filters for p 6 {0, 1, 2} 
and d 6 {1, 2, 3}, while Figure 3.3 shows wavelets 0(x) = ' j*d f° r the same values 
of p and d. Wavelets from this family have a support of length d+p+1, regularity 
order p, and are either symmetric or antisymmetric. 



Table 3.1: Impulse responses h(n), g(n), l(n), and k(n 



n 


h(n) 


d= 1 


d = 2 


d = 3 


g(n) 


l(n) 


k(n) 


9(n) 


l(n) 


k(n) 


g(») 


l(n) 


k(n) 


_2 
















i 






-1 


0.5 


1 






1 






-3 


-0.125 







0.5 


-1 


0.5 


-0.25 


-2 


0.5 


-0.25 


3 


0.625 


0.0625 


1 






0.5 


0.25 


1 


0.5 




-1 


0.625 


-0.0625 


2 


















-0.125 





for p = and d € {1,2, 3}. 



Table 3.2: Impulse responses h(n), g{n), /(??), and k(n) for p=l and d £ {1,2,3}. 



n 


h(n) 
l(n) 


d = 1 


d = 2 


i/(») 


k(n) 


9(n) 


k(n) 


-1 


0.25 


1 


-0.0625 


1 


-0.0625 





0.5 


-1 


-0.3125 


-2 


-0.375 


1 


0.25 




0.3125 


1 


-0.0625 


2 






0.0625 







n 


J = 3 


fe(n) 


g(n) 


l(n) 


fe(n) 


-3 






-0.015625 




-2 




1 


-0.09375 


0.00390625 


-1 


0.25 


-3 


0.265625 


0.04296875 





0.5 


3 


0.6875 


0.1015625 


1 


0.25 


-1 


0.265625 


-0.1015625 


2 






-0.09375 


-0.04296875 


3 






-0.015625 


-0.00390625 



As already discussed, wavelets with p = 2 and d=\ from a family of wavelets 
with p even and d=\ were used in [27, 28], whereas filters with p= 1 and d = 2 from a 
family of filters with p odd and d = 2 were employed by Laine and collaborators [1 1 , 22, 
23]. Here described transform puts no restrictions on the choice of p or d whatsoever, 
and uses a better initialization procedure than the one originally proposed in [28]. 



2(i 



Table 3.3: Impulse responses h(n), g(n), /(»), and k(n) for p — 2 and d (E {1,2,3} 



7? 


fc(n) 


d = 


1 


d = 


2 


9(n) 


l(n) 


k(n) 


9(n) 


'(») 


fc(n) 


-2 


0.125 






-0.015625 






-0.015625 


-1 


0.375 


1 


0.125 


-0.109375 


1 


0.125 


-0.125 





0.375 


-1 


0.375 


-0.34375 


-2 


0.375 


-0.46875 


1 


0.125 




0.375 


0.34375 


1 


0.375 


-0.125 


2 






0.125 


0.109375 




0.125 


-0.015625 


3 








0.015625 









n 


d = 3 


&(n) 


g{n) 


l(n) 


fc(n) 


-4 






-0.001953125 


0.000244140625 


-3 






-0.017578125 


0.003662109375 


-2 


0.125 


1 


-0.0703125 


0.0263671875 


-1 


0.375 


-3 


0.085937 


0.0908203125 





0.375 


3 


0.50390625 


0.13037109375 


1 


0.125 


-1 


0.50390625 


-0.13037109375 


2 






0.0859375 


-0.0908203125 


3 






-0.0703125 


-0.0263671875 


4 






-0.017578125 


-0.003662109375 


5 






-0.001953125 


-0.000244140625 



27 



-0.2 



-0 6 







---71- 


— , 










- 


■ 




]/ ~~\ 


\\ | 










' 


- 


/ 


1 ' 








. 




- 






<' 












\\ 




b 
























-3 


-2 


-1 











2 


3 



(a) 



1 


- 


--jr- 




,__ . 


- 














05 


A 


\\ 




\l 


V 





J 


\\ 




i 


V 






0.5 










- 


-1 






V 7 
\\ // 




- 


1 5 










- 


-2 










- 




Figure 3.3: (a) Wavelets (a) 0(x) = ^ M , (b) 0(ar) = ^ff^ , and (c) 0(*) 
d3/? g; 3 3(j) for p = (dashdotted), p=l (solid), and p = 2 (dashed). 



29 

the zero-crossings of W m s(x) for d = 2, a scale-space image of the representation sim- 
ilar to [48] can be obtained. The only differences are that curves in the scale-space 
plane are computed for dyadic scales only and that i) m is a close approximation of a 
Gaussian instead of a true Gaussian. 



CHAPTER 4 
TRANSFORMS IN TWO DIMENSIONS 

In this chapter, we will extend the one-dimensional transform from Chapter 
3 to two dimensions 1 and derive several two-dimensional transforms with translation 
and at least approximately rotation-invariant decomposition. 

4.1 2-D Discrete Dyadic Wavelet Transform Revisited 

The dyadic wavelet transform of a function s(x,y) G L 2 (R 2 ) is defined as a 
set of functions [28] 

{Wls(x,y),Wls(x,y)} m€Z , (4.1) 

where W^s(x,y) = s * ^ m (x,y) for i = {1,2} and ^ m (x,y) = 2- 2m t l {2- m x,2- m y) 
are wavelets ^'(x,y) expanded by a dilation parameter 2 m . 

To ensure coverage of the frequency space there must exist an A 2 > and 
B 2 < oc such that 

A 2 < £) j:iV'*'(2 m ^,2 m u; v )| 2 < B 2 

m = — oo i=\ 

is satisfied almost everywhere. If (nonunique) functions \ l {x.y). \ 2 (x,y) are chosen 
such that their Fourier transforms satisfy 



£ ^^(2^,2^)^(2^,2-^) = 1. 

l = — oo z'=l 

the function s(x,y) may be reconstructed from its dyadic wavelet transform by 

4^y)= f; s:wi«*xi.(*,y). 

7n=— oo i=l 

where \j n (.r,y) = 2~ 2m x i (2~ m x,2- m y). 

1 For extensions to higher dimensions, please refer to [18, 19]. 

30 



31 

However, when processing discrete functions the scale 2 m may no longer vary 
over all m € Z . Let the finest scale be normalized to 1 and the coarsest scale set to 
be 2 M . Let us, similar to [28], introduce a real smoothing function (f>(x,y) such that 
its Fourier transform satisfies 

|#c*,<*)| a = £ E^(2 m ^,,2"H) \'C2 m ^2 m ^ y ). (4.2) 

m =0 /=1 

Here, as in one dimension, a finite energy discrete function (s(n x ,n y ) £ l 2 (Z )) can 
be written as the uniform sampling of some function smoothed at scale 1: s(n x ,n y ) 
= S f(n x ,n y ), where f(x,y) G L 2 {R 2 ) is not unique, and S m f(x,y) = /* <f> m (x,y). 
This led Mallat and Zhong [28] to a two-dimensional analog of the one-dimensional 
definition of the discrete dyadic wavelet transform: 2 

{ S M -if{n x +s, n y +s),{W^f(n x +s,n y +s),W^f(n x +s,ny+s)} me [o t M-i]} {nxny)&Z 2- 

We will use, as in Section 3.1, a spline-based initialization procedure. 

To implement a multidimensional discrete dyadic wavelet transform within a 
fast hierarchical digital filtering scheme, the wavelets were chosen to be separable 
products of one-dimensional functions: 

1 (*,y) = ^(ar)^(y), (4.3) 

4>\x,y) = V(y) </>(x), (4.4) 

where d(.r) and 0(.r) were chosen as described in Section 3.1 (i.e.. d(x) = /3 p (x) and 
^>(u>) specified by (3.6)). 

From (4.3), (4.4), and (3.6), we may write 



frfatUy) = G- s {u! r )l3 p {u T )l] p {u;y), (4.5). 

j> 2 (u x ,uj y ) = G- S {u y )f3 p {u x )f3 p {uj y ), (4.6) 



2 As in Section 3.1, we put the finest scale of the transform at »7) = 0. 



32 



where G(w) is given by (3.15) for d G {1,2}. Choosing 

A A 

X (u x ,u> y ) = K a (u x )Ti(wy)P P (u) a .)P P (u) y ), 

X 2 (U X ,U}y) = K a (LJy)Tl(U X )P p (u) X )0 p (u}y), 



(4-7) 
(4.8) 



where K(u>) and Ti(uo) are digital filter frequency responses, we may compute (4.2) 
for the finest two scales by 

X> f (2-v,2^) x\2u t ,2u y ) = \Uu x ,ui v )\ 2 - |#(2w„2w,)| 3 . (4.9) 

, = i 

Inserting the terms defined by (4.5), (4.6), (4.7), (4.8), (3.5), and (3.9) with <j>{ut x ,w v ) = 
^ p {u} x )^ p (ijJy) into (4.9) results in 

K{uj x )G{ui x )T x {u v ) + KfaWMTrfe) + \H(u x )\ 2 \H(u y )\ 2 = 1. (4.10) 

Equation (4.10) represents a relation between the frequency responses of the digital 
filters used to implement a multidimensional discrete dyadic wavelet transform and 
is a multidimensional analog to (3.10). 

Solving (4.10) for 2i(u;) by substituting K(uj)G(u)) from (3.10) yields the 
closed formula [28] 

T^)= [ -(l + \H(^\ 2 ) (4.11) 

In Table 4.1 we provide the filter coefficients for Ti(uj) from (4.11) for p £ {0, 1,2}. 
All other filters from (4.10) were already specified in Section 3.1. 

Table 4.1 



: Im 


Dulse responses t 


i(n) for p e 


n 


p=0 


p=l 


p=2 


-3 






0.0078125 


-2 




0.03125 


0.046875 


-1 


0.125 


0.125 


0.1171875 





0.75 


0.6875 


0.65625 


1 


0.125 


0.125 


0.1171875 


2 




0.03125 


0.046875 


•J 






0.0078125 



33 

As in the one-dimensional case, a two-dimensional discrete dyadic wavelet 
transform can be implemented as a fast hierarchical filtering scheme. To derive 
such an implementation, we, similar to the one-dimensional case from Section 3.1, 
use the definition of the two-dimensional dyadic wavelet transform (4.1) and require 
s(uj x ,u} y ) = for \uj x \ > 7r or \u> y \ > k. Using Shannon's sampling theorem in two 
dimensions [16], (4.5), (4.6), and m = 0, we get 



oo oo 



/oo /-co -Q -"-' 

J2 Yl s(i x ,iy)smc(t x -i x )smc(t y -i y y 
-CO J — CO ■ 

l x =-rx., l y = -,yz> 

X 

• Y, 9-s{m)(3 p (x -t x - m)f3 p {y - t y ) dt x dt y 

n 

and 

/CO /-oo ' — ' -'-' 

Y. H s{i x ,i y )smc(t x -i x )smc(ty 
■° oJ -° i x =~ooi v =- o 

oo 

•Pp(x - t x ) J2 9-s(m)P P (y -t y - m) dt x dt y . 



m — — co 



oo oo 



111 = — CO 



As in one dimension, we make use of the fact that the cardinal spline functions 
converge to the sine function as their order r tends to infinity, and write 

W^s(uj x ,u y ) ~ S(u x1 u y )B~ 1 (u} x )B~ 1 (u}y)^ p+r+1 {u} x )^ p+T+1 (u; y )G- a (ut x ) 



an( 



W$s(u} x ,u y ) ~ S(u t ,u y )B r 1 (u x )B r 1 (u}y)$ p+r+ i(u x )^ p+r+1 (u} y )G- s (u; y ), 

or for m £ [0, M) and discrete signal processing 

F{W l m s(x,y) } ~ S(u x ,ujy)B;\u; x )B^(u;y)B p+r+1 (u x y 

x=n x ,y=n v 

m—1 

■B p+r+1 (u y )G^(2 m uj x ) J] H_ s (?u: x )H_ s (?uJ y ). (4.12) 

t'=0 

and 

F{W 2 m s{x,y) _ _ }~S(u x ,u y )B;\u x )B; 1 (u y )B p+r+1 (u, x y 

771 — 1 

■B p+r+1 (u y )G. s (2 m i, y )l[H^C2 i ^.)H_ s C2^ y ). (4.13) 



i=0 



Equations (4.12) and (4.13) describe the decomposition part of the filter bank 
implementation of a two-dimensional discrete dyadic wavelet transform. The recon- 
struction part can be obtained from (4.7)-(4.9) with (3.5) and (3.9). The entire filter 
bank implementation of the transform is shown in Figure 4.1. Except for the pre- 
filtering and postfiltering, we readily recognize the implementation proposed in [28]. 

Using the fact that a wavelet V'( x ) ^ s equal to a first (d= 1) or a second {d = 2) 
derivative of a spline function /? p+ d(x), (4.3) and (4.4) may be rewritten as 

V[x,y) = — -Q^i — , i,d€{l,2], 

where 

i) i (x.y) = (3 p+d (x)My) 

and 

d 2 (x,y) = P P (x)p p+d (y). 

Let us denote W m s(x,y) = (W^s(x, y), W^s(x. y)). V = (£, £), A = V 2 = 
(|^ + Jt-), and assume that d'(x,y, ) can be approximated by O(.v.y) for both / £ 
{1.2}. 

For c/ = 1 it then follows that 

W m s(x,y) = 2 m V{s*# m ){x,y). (4.14) 

Thus for c/ = 2 we can write 

2 



'£Wi nS (x,y) = 2 2m A(s*# m )(x,y). (4.15) 



1=1 



With ${x, y) being a Gaussian, finding local extrema of (4.14) in the direction 
of gradient V corresponds to the filtering stage of a Canny edge detector [6], and 
finding zero-crossings of (4.15) corresponds to the filtering carried out with a Marr- 
Hildreth edge detector (Laplacian of Gaussian) [29]. (Note that both edge detectors 



Bftojlgffly) B p+r+1 (m x ) B^ !(cp y ) 



G-s(co x ) 



G s (co v ) 



H s (co x ) H^) 



(a) 



G,(2co x ) 



G s (2 My ) 



35 



H s (2co x )H s (2o) y ) 



K s (co x ) T^coJ 



1,(0),) K s (o) ) 



Ks(2co x ) 1,(2(0,) 



T 1 (2co x )K s (2 My ) 



'H 



^(2^)^(2^) 



H -sK) H -*(co y ) 



(b) 



B p 1 +r+l (o) x )B p 1 +r fl (cQ y ) -- B r K) B r (co y ) 



Figure 4.1: Filter bank implementation of a two-dimensional discrete dyadic wavelet 
transform (a) decomposition and (b) reconstruction for two levels of analysis. H* s (u) 
denotes the complex conjugate of H- s (uj). 



36 

involve postprocessing). Edge detection based on finding local extrema of W m s(x.y) 
or zero-crossings of £2, =1 W^ n s(x,y) is therefore an approximation to the Canny or 
the Marr-Hildreth edge detector over a range of dyadic scales. The differences stem 
from the fact that r)(x.y) is neither a Gaussian nor is i?*(x,y) equal to i)(x,y). 

4.2 Steerable Functions 

When extending the transform from Chapter 3 to two dimensions, one of the 
first questions that conic to mind is how to deal with the fact that derivatives can 
be defined in any direction of the plane. 3 In case of a first derivative of a Gaussian, 
one can simply compute first derivatives of a Gaussian in x and y directions and 
then determine the derivative in any direction from these two directional derivatives 
[6]. For higher order derivatives of a Gaussian, Freeman and Adelson [12] showed 
that order+1 directional operators are needed for synthesizing the operator at any 
orientation. In fact, functions with the property of expressing their arbitrary rotations 
as linear combinations of some functions are not limited to derivatives of a Gaussian. 
Below, we briefly restate some of the results from [12]. 

A two-dimensional function is called "steerable" if its rotations generate a 
finite dimensional space. Rotations of a steerable function f(i\0) can therefore be 
expressed as 

/M-0 o ) = £ci(0o)/,(r,0), (4.16) 

where 6 denotes an arbitrary angle, {cj(0 o )} stands for a set of interpolating func- 
tions, {/,(r, 0)} is a set of basis functions, and r = y/x 2 + y 2 and 6 = arg(.r,y) are 
polar radius and angle, respectively. 

If /(r, 6) represents a filter kernel, the result of filtering with a rotated fil- 
ter f{r\6 - 6 ) can be computed simply by {c;(0 o )} weighted linear combination of 
outputs from basis filters {f t (r,8)}. Only the outputs from basis filters need to be 



3 Second derivatives of central B-splines can be used, as we saw in Section 4.1, to approximate 
Laplacian of a Gaussian for approximately rotation-invariant, although not directional, processing. 



M 



precomputed and then the output from a filter rotated by any angle 9q can be found 
by interpolating between them. When a large number of rotations of a template filter 
is required, the above scheme can lead to substantial savings in both computational 
cost (time) and memory requirements (space). 

The necessary condition for a function j(i\9) to be steerable is that f(r,9) is 
bandlimited in its polar angle: 

N 



f(r.9) = £ a n (r)e> n6 

n=-N 



it.17) 



This can be verified by inserting (4.17) into (4.16) and by assuming, for con- 
venience, that fi{r,9) = f(r,9 — 9,): 

i 



a n (r)e- jn6 ° =Y / c i (0o)a n {r)e- jnBi , 



(4.18) 



i=i 



where {#,-} is a set of user defined angles and n £ [— A", 0]. 4 Since only nonzero 
coefficients a n (r) are of interest, both sides of (4.18) can be divided by a n (r). This 
yields a matrix equation from which interpolating functions c 4 (# ) can be determined: 



J6o 



.3N60 



=Jd 



JNt 



cJV2 



o]N6 2 



■ J" i 



JNBj 



ci(0o) 
c 2 (0o) 

ci(9 ) 



(4.19) 



For coefficients a n = the rows corresponding to each n were removed from 
the matrix formulation shown in (4.19). For this system to have a solution, the angles 
{9;} must be chosen such that the columns of the matrix are linearly independent. 

In [12] they proved that the minimum number of basis functions /,(;'. 9) needed 
to steer /(r, 9) according to (4.16) is equal to the number of nonzero coefficients a n (r) 
in the Fourier series expansion (4.17). 

To conclude this brief description of steerability, let us only remark that func- 
tions which are not steerable (i.e., do not have a finite number of terms in (4.17)) 



4 Note that the constraints are the same for n £ [-N, -1] and n£[l, .V], so that a subset of all 
possible values for n € [-N, N] can be taken. 



38 

can be approximated with steerable functions (a singular value decomposition was 
employed for approximating such functions efficiently [31]), and that the technique 
of expressing transformed versions of a function as linear combinations of a fixed set 
of basis functions is not limited to rotations (extensions to translations [39], scalings 
[31, 39], and general transformations [14] have been reported). 

4.3 Steerable Dyadic Wavelet Transform 

Returning to the question from the beginning of Section 4.2, the answer seems 
obvious: one needs to construct a steerable analog to the one-dimensional transform 
from Chapter 3. Steerable transforms are nothing new — quite a few [13, 17, 38, 39] 
have been devised, some of them [17, 38] exactly for the computation of directional 
derivatives. Here, we are not interested in any directional derivatives: we want to 
use derivatives of central B-splines which, as the order of B-splines increases, tend to 
derivatives of a Gaussian. 

We define a steerable dyadic wavelet transform of a function s(x, y) G L 2 {R 2 ) 
at a scale 2 m , m 6 Z, as [21] 

Wi a s(x,y) = 8*tj>l(x,y), (4.20) 

where ^(x,y) denotes ^ m (x,y) rotated by 9 U >:\„{x.y) = 2" 2m 0(2-"V,2- m y), 
0(x,y) is a steerable wavelet that can be steered with I basis functions, and t = ^j-k 
with i e {1,2,...,/}. 

Analogously to the one-dimensional case (cf. Section 3.1) we require the two- 
dimensional Fourier plane to be covered by the dyadic dilations of ^ i {2 m u x ,T n u y ): 
there must exist ,4 3 > and B 3 < oo such that 

4» < £ El^'(2 m ^,2 m u;,)| 2 < i?3 (4.21) 

m=— oo i=l 

is satisfied almost everywhere. 



39 

If (nonunique) reconstructing functions \^(:r,y) are chosen such that their 
Fourier transforms satisfy 

f) E ( / > '(2 m u; x ,2 m ^)x ! (2 m ^,2 m u; y ) = l, (4.22) 

m=— oo t=l 

the function s(x, y) may be reconstructed from its steerable dyadic wavelet transform 

by 

oo / 

»{*,V)= E EwL**x\n(*,v)< (4-23) 

m=— oo i=l 

where \ l m {x,y) denotes \ m {x,y) rotated by 6, and \ m (x.y) — 2~ 2m \(2 _m .r, 2~ m y). 

To derive an algorithm for the fast computation of the transform, we, similar 
to (3.3), introduce two smoothing functions such that 

00 I 

i{L> x ,u> y )<p{u x ,u> y ) = £ ^^'(2 m ^,2 m a Jy )i'(2 m u; x ,2 m u; y ). (4.24) 

We choose wavelets that are steerable analogs to the one-dimensional wavelets 
from Section 3.1: 5 

il>(u r ,ue) = (ju} r cos(u e )) d (— s^J , (4-25) 



where u T = «/u;J + ^ and u^ = arg(cj :r ,u; y ). These wavelets can be steered with d+\ 

basis functions (i.e., / in (4.16) is equal to d+1). 

Choosing 

4>(2u r ) - H st (u r )<j)(u) r ), (4.26) 

^ % {u} r ,u 6 ) = G st {uj r ,u: 9 - 9i) 6{u: r ). (4.27) 

ip(2u T ) = L st (u> r )ip(u r ), (4.28) 



and 



X % (w r ,u e ) = K at (L> r ,u> g - Oi)<p{u) r ), (4.29) 



This choice can be viewed as an extension of the transform presented in [20, 21, 24]. 



10 

and using (4.26) through (4.29) with (4.24) computed for the finest two scales, we 
obtain 

^G a t(u; r ,u>4, - O^Kstiur,^ - 0i) + H 3t (ur)L s t{u T ) = 1. (4.30) 

Setting <£(u> r ) = /? p (uv), and employing (3.5) and (4.25) with (4.26) and (4.27), 
we find that 

H«(u r ) = H-fa) (4.31) 

and 

G at {u r ,u e ) = (co S (u e )) d G- 3 {u> r ), (4.32) 

where H(u>) and G(u) are specified by (3.9) and (3.15), respectively. 

By inserting (4.31) and (4.32) into (4.30), the missing two filters can be chosen 
as 

L«(u r ) = L a (u r ) (4.33) 



and 



K at (u r ,U0) = — (cos{u e )) d K s {u r ), (4.34) 



where L{^) and K(u)) are given by (3.16) and (3.17). respectively, and Cd = 
Ef = i(cos(u;-^)) 2d . 

4.4 Multiscale Spline Derivative-Based Transform 

Let us pause here for a brief assessment of the two-dimensional steerable trans- 
form derived so far. We have chosen steerable wavelets (4.25) which are equal to d-th 
order derivatives of circularly symmetric spline functions in the direction of x-axis 
(note that knots for these splines are circles) and we have laid a foundation for fil- 
ter bank implementations in (4.30). An obvious shortcoming of this scheme is the 
fact that none of the filter kernels obtained from (4.31) through ( 1.31) is compactly 
supported on the rectangular grid. For implementations using digital filters, one is 
therefore forced to approximate these frequency responses, and by doing so. (4.30) 



II 

may not hold anymore. Filters in filter bank implementations of steerable pyramids 
described in [17, 38, 39], for example, were designed by using various techniques 
for approximating desired frequency responses. None of the reported filter banks 
achieved perfect reconstruction and all filter kernels were nonseparable. Here, we 
will take a different approach. We will approximate the wavelets in (4.25) in a way 
that will lead to x-y separable filters in a perfect reconstruction filter bank imple- 
ment at ion of the transform such that the quality of approximation will improve by 
increasing the order of B-splines. 

Let us approximate wavelets from (4.25) with 

, v d (3 p+ d(x) 

nx,y) = — -jp — i3 P +d(y)- (4.35) 

Based on the fact that B-splines tend to a Gaussian as their order increases, 
it is easy to see that both wavelets (4.25) and (4.35) converge to the same functions 
(i.e., d-ih order derivatives of the normalized Gaussian in the direction of .r-axis) as 
p — * 'CO. 

In order to steer wavelets tj){x,y) given by (4.35) (note that steering will 
be only approximate, since these wavelets are not steerable), we need to find basis 
functions that will approximately steer ij)(x,y). Computing rotations, as we did in 
(4.20), is not practical here, because arbitrary rotations of (4.35) cannot be expressed 
exactly in terms of central B-spline functions from Chapter 2. Instead, we take 
advantage of the property of circularly symmetric functions that rotations of their 
directional derivatives are equal to directional derivatives in rotated directions: 

d d Q c (x,y)\ d d g c (x,y) 



Tin 



dn d dni 



"0 



where 1Z 6q stands for rotation by 6 , dec £' v) = nVg c (x,y), g c (x,y) is a circularly 
symmetric function, and ng denotes vector n = (cos 0, sin 6) rotated by 6 . 
Let us choose 

g(x,y) = 8 p+d (x)/3 p+d {y), 



12 



which is approximately circularly symmetric function for higher order splines. A 
rotation of %l>(x,y) = f^J y from (4.35) by 6q can therefore be approximated by 

where n = (cos $o, sin #o) = (n x ,n y ). 

Note that in case of Gaussian, which is both x-y separable and circularly 
symmetric, (4.36) becomes exact (e.g., for g(x,y) = e~^ x + y ~>, o = —0, and d = 
{2,4}, we obtain, up to a scaling factor, x-y separable basis set for the second and 
fourth derivative of Gaussian from Tables III and VII of [12]). 

Having found a set of basis functions (4.36) that approximately steer wavelets 
(4.35), we want to construct a transform such that Equations (4.20) through (4.24) 
will be valid (superscript i must be viewed now as an index, rather than rotation by 
0i). In frequency domain, we can express basis functions from (4.36) as 

4> i+1 (u x ,u; y ) = G±?(u x )GL,(Ljy)P P+i {u x )P P+i - i (u y ), £ = 0,1 d, (4.37) 

where G d (u) is given by (3.15) and G°(w) = 1. 

Choosing appropriate x % (u x ,u> y ) to obtain a relation needed for the filter bank 
implementation of the transform is more complicated than in one dimension. Since we 
could not find a general solution for arbitrary d, we solve each case separately. Below, 
we present solutions for the first four orders. When d < 2, we impose ^(w^u,,) = 
<f>(u> x ,u}y) = /3 p (u> x )fl p (u}y), a constraint analogous to the one from Section 3.1. 

For d = 1 , we write similar to [28] 

x\u>*,Uy) = K]{u x )T l {uj y )i3 p (uj x )fip-\{uy). (4.38) 

\ 2 (u x ,^ y ) = r 1 (w x )j£(<«;J&_ 1 (w x )&(a; y ), (4.39) 

where K d (u>) and Ti(u) are given by (3.17) and (4.11), respectively. 



13 

Computing (4.24) for the finest two scales and inserting (3.5). (4.37), (4.38), 
and (4.39) yields a relation between frequency responses 

G\u x )K\u> x )Ti{u>y) + TiK^Vv^Vv) + \H(u t )H{uy)\ 2 = 1. 
For d = 2, we choose 

\\u; r ,uj y ) = K*(u x )T 2 (uJy)Pp(u x )P P -2(u y ), (4.40) 

X 2 {u x ,u v ) = K](u x )K}((jj y )P P - 1 (u g )P P -.i(u y ), (4.41) 

X 3 (w*,w y ) = T 2 (u; ;r )A' s 2 (u;j / )^ p _2(u; ;c )^ p (LJ y ), (4.42) 



w 



here 



r a (w) = |#(u,)| 2 . (4.43) 

Using (3.5), (4.37). and (4.40) through (4.42) with (4.24) results in 

G 2 (u x )K 2 (u x )T 2 (u y ) + G\u; x )K\u x )G l {w y )K\u y ) + T 2 (u x )G 2 (u y )K 2 (u y )+ 

+\H{u x )H(u> y )\ 2 = l. 

For orders d > 2, we require <f>(u> x ,u y ) = J3 p (u} x )P P {u> y ) and ip(u x ,u> y ) = (f>(u} x )<p(uj y ), 
where fi{u>) is specified by (3.7) and (3.16). 

For d = 3. we choose reconstructing functions 

^(u^Uy) = K*(w x )<p(u x )<p- 3 (L,j y ), (4.44) 

X 2 (w„w y ) = -A' s 2 (^)A' s 1 (cj y )V3(u; r )V3(^)^-i(^ x )^-2(^), (4-45) 

fK,w s ) = -A s 1 (^)A 2 (u; y )V3(a J;r )l / 3K)^-2(^)^-iK). (4.46) 

X*(ws,w y ) = K*(u y )<p- 3 (u x )<f(u) y ), ( 1.47) 

where 

7 3 (w) = -^(1 - \H(lo)\ 2 ), (4.48) 

and £_,-(a>) € L l (R) denotes a function such that <p(uj) = /3,_i(u;)<£_ t -(a;), i € iV. 



II 

Employing (4.37), (3.5), (3.7), and (4.44) through (4.47) with (4.24) yields a 
relation 

-G\u, x )K\u; x )V 3 (u; x )G\u: y )K 2 (u y )V 3 (u Jy ) + G 3 MK 3 M+ 

+H(u x )L(u> x )H(u y )L(u y ) = 1, 

where L(*c) is specified by (3.16). 
For d = 4, our choices are 

X l (e>z,w y ) = Kt(oj x )T 2 {uj y )^(u) x )^- 4 M, (4.49) 

X 2 (wx,u y ) = K 3 (u x )Kl(u y )<p-i(u x )<p- 3 (u> y ), (4.50) 

X 3 (w*» w i») = - K Uvx)Kl(uy)V 4 (u x )V4(u y )<p-2(ux)<P-2(vv), (4.51) 

X 4 (u x ,u y ) = A' s 1 (^ x )A' s 3 (u.' y )^_ 3 (u; J .)^_i(^ !/ ), (4.52) 

X B (w t ,Uy) = T 2 (u> x )K 4 (u> y )<p- 4 (u x )ip(u v ), (4.53) 



where 



V 4 (u x ) = 1 - |tf(u;)| 2 . (4.54) 



Using the above functions (4.49) through (4.53), (4.37), (3.5), and (3.7) in 
(4.24) computed for the finest two scales gives 

G 4 (u x )K 4 (u x )T 2 (u y ) + G 3 (u; a .)A' 3 (-.v)G 1 (u; J ,)A' 1 (u; J ,)- 
-G 2 {u x )K 2 (u x )V 4 (u; x )G 2 {u y )K 2 {u; y )V 4 (u y ) + G\u^ 

+T 2 {uj x )G 4 {u; y )K 4 {uj y ) + H(u x )L{u x )H(u y )L(u> y ) = 1. 

Here, we have even more freedom for choosing the reconstructing functions 
than in one dimension. The above functions for d = {2,3,4} were found by trying 
to imitate the one-dimensional transform from Chapter 3 as much as possible. All 
decomposition filters G d (m) were first paired with corresponding reconstruction filters 



45 

K d (uj), and then all other potential digital filters were specified as polynomials in 
H- a (u). We inserted thus specified filters into a relation between their frequency 
responses and solved for the unknown polynomial coefficients. Since we allowed more 
filters with higher-degree polynomials than expected in the solution, the resulting 
system of linear equations was underdetermined. This allowed enough freedom for 
removal of undesired digital filters and for balance between degrees of polynomials. 

The described procedure for determination of reconstructing functions and 
filters involves quite a lot of heuristics to obtain the appropriate solution from the 
underdetermined linear system. Unfortunately, we are not aware of any systematic 
way (aside from numerical optimization, which may be pretty cumbersome) to obtain 
solutions comparable to the ones above. 

Next, we will derive a filter bank implementation of the transform. As in 
Section 4.1, we assume a bandlimited input signal: s(uj x ,u} y ) = for |o; z | > it or 
|a> y | > 7r. Using Shannon's sampling theorem in two dimensions [16] with (4.20) and 
basis functions from (4.37), we can write 



/OO /"CO ^ -~ 

J2 J2 8 (**i »v)mhc(** - i x )ainc(t y - i y )- 
-OO J — OO ■ 

l x = — OO ly = — OO 

OO OO 

Y. g d -7( m ?){3p+i(x -t x -m x ) Y. g- s {m y )f3 p+d -i{y -ty-m y )dt x dt y , 



m x = — oo m,u=—oo 



where i = 0, 1, ...,(/ as in (4.37). 

Again, we approximate sine functions with r-order cardinal splines, then use 
(2.4), and get 

x — n x ,y — fly 

m—\ 

■B p+r+d . i+1 (u; y )G d - i (2 m u; x )GU2 m u y ) J] H p _ + a i (2 n u x )H^ d - i (2 n u; y ). (4.55) 

71=0 

Using (4.55) with an approximation 5 p+r+t+1 (w) ~ B p+r (uj)Bi(uj), we can ob- 
tain a filter bank implementation of the transform decomposition. The reconstruction 
part follows from (4.24), (4.37), and reconstructing functions for distinct values of d. 



46 

Figure 4.2 shows filter bank implementations of a multiscale spline derivative-based 
transform for d — {1,2,3,4}. For d = 1, we recognize (except for the prefiltering 
and postnltering) the filter bank implementation of a two-dimensional discrete dyadic 
wavelet transform from [28]. For d = 2, however, our transform differs from the filter 
bank presented in [22] (i.e., the corresponding transform described in Section 4.1): 
second derivative is computed only in the directions of x and y-axis in [18, 22], which 
is not enough for steering. Although not particularly appropriate for orientation anal- 
ysis, such a scheme may still, as we have seen in Section 4.1, efficiently approximate 
Laplacian of Gaussian across dyadic scales. 

A transform similar to the one described in this section, was presented in 
[17, 38, 39]. Their filter bank implementation is less redundant (downsampling is used 
on the lowpass branch, while simultaneously keeping aliasing negligible by using a 
filter with insignificant amount of energy for \uj r \ > |) and uses reconstruction filters 
with same magnitude frequency responses as the decomposition ones — a possible 
advantage for certain applications. They have, on the other hand, little control over 
the function from which derivatives are computed (to obtain a d-th derivative, they 
multiply a circularly symmetric filter by (-jcos6) d with all filters being obtained 
by recursive minimization of a weighted combination of constraints), filter bank does 
not have perfect reconstruction, and none of the filters is x-y separable. 



IT 



B> x ) B> y ) -JJ^hK) B p+r+1 (co y ) 

(a) 



B p+r+1 (ro x ) B p+r+1 (co y ) B r (o) x ) B r (co y ) 
(b) 



G 1 s (2 m co x ) 



G 1 s (2 m co v ) 



, sV - w y , 



H s (2 m (o x )H,(2V) 



(c) 



G 2 (2 m co x ) 



Qs(2 m a) x ) G 1 s (2 m co y ) 



G 2 (2 m co Y ) 



H,(2 m co x )H,(2V) 



(e) 



Kj2 m a> x ) T,(2 m co v ) 



"1 \ TT-l/^m 



1,(2©,) K;(2 m co v ) 



L s (2 m o) x ) L s (2 m co y ) 



(d) 



K 2 (2 m (o x )T 2 (2 m co v ) 



•'x/ ^V* - "'y/ 



K;(2 m co x )K;(2 m o) y ) 



T 2 (2 m co x ) K 2 (2 m co v ) 



L s (2 m co x ) L s (2 m co y ) 



(f) 



Figure 4.2: Filter bank implementation of a multiscale spline derivative-based trans- 
form for m G [0, M - 1]: (a) Prefiltering, (b) postfiltering, (c) decomposition and 
(d) reconstruction modules for d = 1, and (e) decomposition and (f) reconstruction 
modules for d = 2. 



is 



G 3 s (2 m co x ) B^V) 



G 2 (2 m (o x )G.;(2 m (o y ) 



Gl(2 m co x ) G 2 (2 m co Y ) 



B 2 (2 m co x )Q 3 s (2 m (o y )- 



H s (2 m (o x )H,(2 m co y ) 

(g) 



K^coJ EfcV) 



-K 2 (2 m (o x )V 3 (2 m co x )K;(2 m co,)V3(2 m (o y ) 



K;(2 m (o x )V 3 (2 m (o x )K 2 (2 ni co y )V3(2 m co y ) 






^(2 m co x )K 3 (2 m co y ) 



L s (2 m co x ) L S (2V) 



(h) 



G 4 (2 m (0 x ) B,(2 m co y ) 



G^X^coJB^V) 



G 2 (2 m co x ) G 2 (2 m (o v ) 



G^coJ B 2 (2 m a) x ) C&V) 



B 3 (2 m (o x ) G 4 (2 m co y ) 



H :S (2 m o) x )H s (23) 



(i) 



K 4 (2 m (o x )B 3 , (2 m (o y )T 2 (2V) 



K 3 (2X)B,\2V)K:(2 m co y ) 



K 2 (2 m (o x )V 4 (2 m (o x )K 2 (2 m (o y )V 4 (2 m (o y ) 



R 2 , (2 m a) x )K;(2 m co x )K 3 (2 m (o y ) 



B 3 1 (2 m (o x )T 2 (2 m (o x )K 4 (2 m (0 y ) 



1 LsPX) L s (2 m co y ) 
0) 



Figure 4.2: Continued: (g) Decomposition and (h) reconstruction modules for d = 3, 
and (i) decomposition and (j) reconstruction modules for d = 4. Decomposition 
modules are recursively applied at the locations of the filled circles. 



CHAPTER 5 
IMPLEMENTATION ISSUES 

Except for the steerable dyadic wavelet transform, all transforms presented in 
Chapters 3 and 4 can be implemented as filter banks consisting of one-dimensional fil- 
ters only. In this chapter, we show how to take advantage of symmetry/antisymmetry 
of filters when combined with a mirror extended input signal. 

5.1 Finite Impulse Response Filters 

Since all two-dimensional filters used in the filter bank implementations of 
the transforms from Chapter 4 are x-y separable, we will begin this section with a 
detailed description of FIR filter implementations for the one-dimensional discrete 
dyadic wavelet transform implementation from Figure 3.1. The extension to two 
dimensions will then be straightforward. 

Let us refer to filters applied at scale 2 m as filters at level m+1, and let filters 
at level 1 (Equations (3.9), (3.15) through (3.17), (4.11), (4.43), (4.48), and (4.54)) 
be called "original filters,'' to distinguish them from their upsampled versions. As 
an input to the filter bank from Figure 3.1, we consider a real signal s(n) G 1 2 {Z). 
n G [0, N - 1]. Depending on the length of each filter impulse response, filtering an 
input signal may be computed either by multiplying the discrete Fourier transforms 
of the two sequences or by circularly convolving s(n) with a filter's impulse response. 1 
Of course, such a periodically extended signal may change abruptly at the boundaries 
causing artifacts. A common remedy for such a problem is realized by constructing 
x As is customary in image processing, we use circular, rather than linear, convolution. 



1!) 



50 
a mirror extended signal [18] 

s (n) = l S{ - n - l) if « € [-iV, -1] 
mel j \ s(n) line [0,7V -1], ( ° j 

where we chose the signal s me (n) to be supported in [— TV, N — 1]. It will become 
evident shortly, that mirror extension is particularly elegant in conjunction with 
symmetric/antisymmetric filters. 

Let us first classify symmetric/antisymmetric real even-length signals into four 

types [30]: 

Type I f(n) = /(-»), 

Type II /(n) = f(-n - 1), 

Type III /(n) = -/(-«), 

Type IV /( n ) = -/(-n-l), 

where n € [-A r , A' - 1]. Note that for Type I signals the values at /(0) and f(—N) 
are unique, and that for Type III signals the values at /(0) and f( — N) are equal to 
zero. 

Using properties of the Fourier transform it is easy to show that the convolu- 
tion of symmetric/antisymmetric real signals results in a symmetric/antisymmetric 
real signal. If a symmetric/antisymmetric real signal has an even length, then there 
always exists an integer shift such that the shifted signal belongs to one of the above 
types. 

Now, we are ready to examine the filter bank implementation of a one-dimensional 
discrete dyadic wavelet transform from Figure 3.1 with filters given by (3.9) and (3.15) 
through (3.17) driven by a mirrored signal s me (n) at the input. Let the number of 
levels M be restricted by 

M < l + log 2 / ~* (5.2) 

■'-'max -1- 



51 

where L max is the length of the longest original FIR filter impulse response. 

Each FIR filter block in the filter bank consists of a filter and a circular shift 
operator. Equation (5.2) guarantees that the length of the filter impulse response 
does not exceed the length of the signal at any block. 

Since our input signal s me (n) is of Type II and noninteger shifts at level 1 are 
rounded to the nearest integer, it follows that a processed signal at any point in the 
filter bank belongs to one of the types defined above. This means that filtering a 
signal of length 2 A* can be reduced to filtering a signal of approximately one half of 
its length. (For Types I and III, A + 1 samples are needed. However, for Type III 
one needs to store only A — 1 values because zero values are always present at the 
zeroth and (— JV)-th sample position). 

Implementation is particularly simple for FIR filters designed with d even and 
p odd. Filters are of Type I in this case, so the signal at any point of the filter bank 
will be of Type II. An FIR filter block from the filter bank shown in Figure 3.1 can 
therefore be implemented by 

F s , m u(n) = f(0)u II (n) + Y,f(i)[u II (n-2 m i) + u JI (n + 2 m t)l n € [0, A r -1], (5.3) 
where 



u//(n) = < 



u(-n-l) ifne[-f,-l] 

u(n) ifn€[0,AT-l] (5.1) 

u(2N-n-l) ifneJA,^-], 



u(n) is an input signal to a block, /(??) is an impulse response of some original filter. 
L is the length of the filter, and A* is the length of an input signal s(n) to the filter 
bank. Implementation of filters b p (n) used for prefiltering and postfiltering represents 
a special case of (5.3) with m = 0. 

A filter bank with the above implementation of blocks and signal s(n) at the 
input yields equivalent results as circular convolution for s me (n) as defined by (5.1). 
In addition to requiring one half the amount of memory, the computational savings 



52 

over a circular convolution implementation of blocks are, depending on the original 
filter length, three to four times fewer multiplications and one half as many additions. 

A similar approach can be used for other filters. However, things get slightly 
more complicated in this case, because the filters are not of the same type and the 
signal components within the filter bank are of distinct types. As a consequence, an 
implementation of blocks that use distinct original filters may not be the same, and 
the implementation of blocks at level 1 may differ from the implementation of blocks 
at other levels of analysis. 

The decomposition blocks at level 1 can be implemented by 

2 1 

G- a , u(n) = J2 g{i)[u n {n - i - 1) - u n (n + i)], n G [1, JV - 1], 
for d odd, (5.3) for d even, 

H- s ,ou{n) = Y, h{i)[un{n-i-l) + u n (n + i)], n G [O.A], 

t'=0 

for p even, and (5.3) for p odd, where «//(/) is defined by (5.4), g(n) and h(n) are 
impulse responses of the filters computed from (3.15) and (3.9), respectively, and L 
is the length of the corresponding impulse response. 

The output from a block G- 3 (u>) at level 1 is of Type III for d odd and of 
Type II for d even, while the output from H_ a (u) at the same level is of Type I for 
p even and of Type II for p odd. 

The decomposition blocks at subsequent levels m G [l,Af — 1] can be imple- 
mented by 

G- a<m u(n) = £ g(i)[ui(n - 2 m (i + s)) - Ul (n + 2 m (i + s))}, n G [l,iV- 1], 

;=0 

for d odd and p even. 

2 l 

G- a , m u(n) = Y, 9(i)[un(n - 2 m {i + s)) - u n (n + 2 m (i + s))}, n G [0, N - 1], 

f=0 



53 



L-\ 
2 



for d and p odd. 

F_ s>m «(n) = /(O)ur(n) + £ /COM" " 2™*') + "/(« + 2 "'0]- n € [0, AT], (5.5) 

1=1 

with f(n) = g(n) for d and p even, 

H.,, m u(n) = J2 h(i)[ Ul (n - 2 m (i + a)) + u/(n + 2 TO (* + a))], n € [0, AT], (5.6) 



i=0 



for p even, and (5.3) for p odd, where 

»(-??) 
u/(n) == { u{n) 



ifn€[-f,-l] 

if n € [0, N\ 



(5.7 



u(2N-n) \ine[N + l,m 



The outputs from blocks G- S (2 m u:) are of Type III for d odd and p even, of 
Type IV for d and p odd, and of Type I for d and p even, whereas the outputs from 
H- S (2 m u)) are of Type I for p even and of Type II for p odd. 

Next, the reconstruction blocks at level 1 can be implemented by 



K s .ou{n) = J2 Hi)[uin(n -?' + !)- u ni (n + i)], n € [0, N - 1], 
i=i 

for d odd, (5.3) for r/ even, 

L 

2 

L,, u(n) = £ '(»)[«/(» - »' + 1) + "/(» + *)]i n G [0, A r - 1], 

! = 1 

for p even, and (5.3) for p odd, where 
uin(n) = < 



-u(-n) 


ifne[-f. 


-11 





if n = 




u(n) 


if n e [1,N 


-1] 





\(n = N 





[5.8) 



( -u(2tf-n) ifne [W+l,^], 



«/(/?) is as defined by (5.7) and k(n) is an impulse response of the filter from (3.17). 
Note that both outputs from blocks K a (u) and L a (u) are of Type II. 

The reconstruction blocks at subsequent levels can be implemented by 

i-l 

2 * 

K a<m u(n) = £ H* + !)[«///(» - 2 m (i + s)) - U//7 ( n + 2 m {i + .s))], n <= [0, AH, 

!=0 



54 

for d odd and p even, (5.5) with f(n) = k{n) for d and p even, 

A' s , m »(n) = £ fc(i + l)[u/v(n - 2 m (i + s))- u IV (n + 2 m (i + s))], n € [0. N - 1], 

i=0 

for d and p odd. 

L sjn u(n) = #_ s , m u(n), 

for p even, and (5.3) for p odd, where «///(/) is given by (5.8), 

' -«(-„- 1) ifne[-f.-l] 

u /v (n) = I n(n) if n € [O.N - 1] 

, -«(2JV-n-l) ifnG[A',f-], 

and H- S _ m u(n) is specified by (5.6). We observe that the outputs from blocks K s ('2 m uj) 
and L s (2 m u;), m G [1. M — 1], are of Type I for p even, and of Type II for p odd. 

When we compare the above implementation of blocks to circular convolution 
driven by a mirrored signal s me (n) at the input, we observe that approximately 
twofold less memory space, three to four times fewer multiplications and one half 
as many additions are required. (For Type I signals an additional sample has to be 
saved because two values are without a pair). 

Until now. we have talked only about the one-dimensional case, whose filter 
bank implementation is depicted in Figure 4.2. Two-dimensional transform filter 
bank implementations (Figures 4.1 and 4.2) are not only comprised of x-y separable 
filters solely, but also use all the filters from Section 3.1. Everything presented in 
this section so far is therefore directly applicable to the two-dimensional case. Filters 
which have nol been t reated yet (i.e.. /,(//). /_>(;/). r : ,,( n ). r,( n ). and filters b p (n ) from 
the decomposition modules of Figure 4.2) can all be realized by (5.3) for p odd or 
m = and (5.5) otherwise (/(??■) denotes an impulse response of any of the above 
mentioned zero-phase filters in this case). 2 



In case of filters i>3(n) and 1^(71), a slightly more efficient implementation can be obtained 
by precomputing new filters k * i>3<n) and k * 1)4(71), and then implementing them by K s m u(n), 
me [0,M- 1]. 



I 1 u- implementation presented in this section performs all operations in the 
spatial domain, however, one could also implement the structures shown in Figures 
3.1, 4.1, and 4.2 with an input signal s me (n) (Equation (5.1)) in the frequency domain. 
For short filter impulse responses, such as those given in Tables 3.1. 3.2 and 3.3, the 
spatial implementation described in this section is certainly more efficient. For long 
filter impulse responses, however, filtering is faster if implemented in the frequency 
domain. Additional details on alternative implementation strategies can be found in 
[33]. 

5.2 Infinite Impulse Response Filters 

Implementation of IIR filters b' 1 ^) which were introduced in Section 2.2 is 
a bit more involved than the one encountered in the previous section. Fortunately, 
the number of different cases is much smaller here: possible input to b^ln) in filter 
banks from Figures 3.1, 4.1, and 4.2 is either of Type II or of Type I. 3 We will use 
ideas and a few results from [43]. 

Let us first take a closer look at the system function B^iz). This function 
can be written as a cascade of terms 

E ^ - j^h^ - ,i- a --.) ( i- a -, - ^ 

which can be expressed in a parallel form as 

-a / 1 1 

-a 2 Vl -or 1 + Y^az 



E ^ = ^-i 7— nr + i-^rr-0- <wo) 



where a and - are poles of the causal and the anticausal filter, respectively. 
The impulse response of this term can be written as 



1 — cr 



3 Note that symmetry types in this section slightly differ from those defined in Section 5.1: here, 
mirror extended signals are periodically repeated, so that they stretch from — oo to oc. 



.-,(, 



We choose to implement E{z) in a cascade form and therefore extract the 
difference equations from (5.9): 



c + (n) = u(n) + ac + (n - 1) n = 1,2 V-l. 



(5.11) 



and 



c(n) = a(c(n + l)-c + (n)) n = N-2,N-3,. . .,0, (5.12) 

where u(n) denotes the input to the block, c + (n) is the output from the causal part, 

and c(n) stands for the output from the block. 

To solve (5.11) and (5.12) we need boundary conditions c + (0) and c(N — l). 

Let us begin with filters fcr 1 (n) in filter bank implementations from Figures 3.1, 4.1, 

4.2(a). 4.2(b). and Figures 4.2(h) and 4.2(j) with m = 0. We derive 

■ ' v_1 a' +l + o 

c+(0) = £ a-u„,(0 = «(0) + £ - T f^ 

:'=— oo i=0 

and, using parallel form (5.10) 



io 



v27V 



-</(/ 



u(0) + 5> <+1 ti(«), (5.13) 



i'=0 



c(N-r 



—a 
1 -a' 2 

— a 

1 -a 2 



( c + (iV-i)+]T 



-i n -V-, + Q .v+i +i 



I 



,2V 



u(iY 



i=0 
N-l 

[c + (N-l)+ Yl ^-^(i)), 

i=.V — 1 — i'o 



(5.14) 



w 



here 



u 1Ip (n) 
»//(") : 



«//(;? mod (27V)) if n>0 

»//(-(n + l) mod (2JV)) if n < 0, 

u(n) if n e [0, A r - 1] 

u(2JV-n- 1) if/? € [.V.2.V- 1], 

iV is the length of an input signal to the filter bank, and ?o < N — l is selected such 
that a la falls below a predefined precision threshold. 

For IIR filters from Figures 4.2(h) and 4.2(j) with m G [1. M - 1] and p odd, 
we get 



o 
c+ (°) = £ [a-**u IIp (i 

t=— oo 

oo N-l 

Y.Y. 

n=0 i=0 



OO iV— 1 .. r -, N 

.„. W--V ^^ I r »+l+2Nn 1 2.Y(n+l)-i | 

= "(0) + EE [o^ 5 ^] + a-^r— u (i) 



(5.15) 



57 



and 



N-\ 



^(A f -i) = T^ T (^(-v-i) + EL{h 



2?n 



+ 



n—0 i'=0 



Q 



.Y+l + i + 2.Yn 



}«(0, (5-16) 



where 



[«1u - 



n 



x if .r € Z 



otherwise. 

If N is a power of two and 2 m_1 < JV ((5.2) guarantees that the latter condition 
is always true) (5.15) becomes 



c+(o) = «(o)+ Y. 



i+1 

a 2 m 



+ 



2-V-i 



Q 2" 



;=o 



W 



hile (5.16) can be written as 



c(X-l) = -^(c + (N-l)+Y 

— o ^~^ 



1 — Q 2™ 



N-i 
Q 2 m 



a «(0, 



+ 



■ v tH' 



Q 2 m 



1 - Q2 m 



-»(/)• 



Finally, the boundary conditions for filters ^(rc) from Figures 4.2(h) and 
4.2(j) with m E [1,M — 1] and p even are 



c+(0)= £ [a-^u/pCi) 

i=— oo 

V^ f r 2ATn 1 

^ j [a -5™-] i/(0)+ a^ 7- 
+ 



n=0 
N-l 



+ 2JVn 



2A r |7i + l)- 
Q 2"> 



u(N)+ 
u(i)\ 



(5.17) 



and 



where 



—a 



c(7V) = - -(2c + (AT)-«(JV)), 



1 -a 



uj(n) = 



«/p(«) = »/(H mod (2A*)). 
' u(n) if n € [0. A"] 



u(2JV-n) if n e [N + 1,2N-1], 
and c~(N) = c + {N) with c~(n) denoting response of the anticausal filter from (5.10) 
was used. 



oS 



Again, if N is a power of two and 2 m l < N, (5.17) can be simplified 



+ m\ — 



c + (0) 



u(0) + ai^u(N) + T.l^{ 



Q2 r 



+ 



2iV-i 



}tl(l 



1 2i£ 
1 — a 2™ 



Series in expressions for c + (0), c(N — 1). and c(N) for filters from Figures 
4.2(h) and 4.2(j) with m G [1. M — 1] can be, similar to (5.13) and (5.14), truncated 
according to the desired precision. 

For orders p greater than three, we implement B~ 1 {z) as a cascade of terms 
E{z) with different q's. Note that the output from block E(z) is always of the same 
type as the input to it. 



CHAPTER 6 
APPLICATIONS 

In this chapter, we present a. comparison between multiscale spline derivative- 
based transform and traditional wavelet techniques on two sample applications. 

6.1 Image Fusion 

Image fusion combines particular aspects of information from the same imag- 
ing modality or from distinct imaging modalities and can be used to improve the 
reliability of a particular computational vision task or to provide a human observer 
with a deeper insight about the nature of observed data. Whether it is combining 
different sensors or extending the dynamic range of a single sensor, the goal is to 
achieve more accurate inferences that can be achieved by a single sensor or a single 
sensor setting. 

The simplest method of fusing images is accomplished by computing their 
average. Such a technique does combine features from input images in the fused 
image, however, the contrast of the original features can be significantly reduced. 
Among more sophisticated methods, multiscale and multiresolution analyses have 
become particularly popular. Different pyramids [5, 41] and wavelet-based techniques 
[7, 20, 26. 32] have been applied to this problem. Fusion is performed in the transform 
space by computing local statistics across the decomposition scales. Typical size of 
neighborhood is between a single pixel and 5 x 5 area with some loss of contrast 
reported for the latter [5]. Some authors argue about of the particular criteria for 
fusion [5, 26], while others propose a set of different combination rules to perform 
image fusion in a variety of ways [7]. 



:»!) 



60 

Here, we are not trying to find the best criteria for fusion. We simply use 
a "maximum magnitude" rule (i.e., at each position and scale of the transform the 
corresponding transform coefficients are compared and those with greater magnitude 
included for reconstruction) and compare the results obtained by traditional wavelet 
techniques with the ones produced by the multiscale spline derivative- based trans- 
form. 

Figure 6. 1(a), (b) shows a pair of images with distinct (lower in Figure 6.1(a) 
and upper in Figure 6.1(b)) areas blurred. Figure 6.1(c) demonstrates the ideal 
result of image fusion obtained by manually cutting and pasting the upper part of 
Figure 6.1(a) and the lower part of Figure 6.1(b). Image fusion resulting from the use 
of discrete wavelet transforms is shown in Figure 6.1(d) for the orthogonal wavelet 
DAUB12 [7] and in Figure 6.1(e) for the linear phase biorthogonal wavelet Bior6.8. 
Figure 6.1(f) illustrates the result of image fusion with a multiscale spline derivative- 
based transform being employed. Note that the image fused using the multiscale 
spline derivative-based transform is virtually indistinguishable from the one obtained 
manually. Both orthogonal and biorthogonal wavelet transform exhibit artifacts. On 
the average, biorthogonal transform performs better than orthogonal due to the linear 
phase of wavelets and filters, but still causes serious artifacts. 

Basing the decision rule on an area rather than on a single pixel certainly 
reduces artifacts of traditional wavelet methods, however, it reduces sharpness as 
well. As noted by Burt and Kolczynski [5], "there is a tradeoff between sharpness 
and shift invariance." Here, we argue the point that the proper choice of transform 
can eliminate artifacts due to translation and rotation invariance altogether. 

Another possible source of artifacts is misregistration. Low-level fusion al- 
gorithms typically require input images that have already been properly aligned. 
Misregistration causes artifacts [7], and it can be recognized from our discussion 



1,1 

in Section 1.1 that they are more obvious for translation and rotation-noninvariant 
wavelet techniques than for the multiscale spline derivative-based transform. 

6.2 Image Enhancement 

Similar to the previous section, we will point out problems associated with 
traditional wavelet analysis on an example. Below, we briefly present the problem of 
contrast enhancement in digital mammography. 

In mammography, early detection of breast cancer relies upon the ability to 
distinguish between malignant and benign mammographic features. Contrast en- 
hancement can make more obvious unseen or barely seen features of a mammogram 
without requiring additional radiation. The complexity of mammographic images 
and the subtlety of malignancies can present a challenge even to expert radiologists. 
In addition to dealing with barely visible mammographic features, human observers 
sometimes simply overlook abnormalities. Such mistakes affect the number of false- 
negative cases considerably. Computer processing of digital mammograms can assist 
radiologists to reach a correct diagnosis more consistently. A variety of computer 
based techniques have been reported in almost three decades of research. The advent 
of direct digital mammography devices has made digital image processing techniques 
more attractive for screening. 

Contrast enhancement in a wavelet transform framework can be achieved by 
applying a (typically nonlinear) function to wavelet coefficients and then reconstruct- 
ing an enhanced image with modified coefficients. We are not trying to find the best 
method for a specific type of images or noise here, rather we want to illustrate pos- 
sible sources of artifacts when orthogonal and biorthogonal wavelet transforms are 
used. 

Figure 6.2(a) shows an original mammographic image. Figure 6.2(b) depicts 
the result of enhancement using piecewise linear enhancement function [11, 23] to 



62 







(c) 



(d) 





Figure 6.1: (a) An image with lower part blurred, (b) An image with upper part 
blurred, (c) Fused image obtained by combining images from (a) and (b) man- 
ually. Results of fusion performed by (d) discrete orthogonal wavelet transform 
(DAUB12), (e) discrete biorthogonal wavelet transform (Bior6.8), and (f) multiscale 
spline derivative-based transform. 



63 

modify wavelet transform coefficients, and Figure 6.2(c) shows the result of enhance- 
ment using the multiscale spline derivative-based transform and the same enhance- 
ment function as Figure 6.2(b). In this figure, we. more than trying to point out 
better enhancement results in Figure 6.2(c). observe that Figure 6.2(b) shows fea- 
tures (artifacts) that are not present in the original image (they are particularly 
obvious in the area right from the mass in Figure 6.2(b)). Here, the goal of enhance- 
ment is making visible features that are present in the original image but obstructed 
by different structures within the same image. Generation of artificial features, while 
not being rare when nonlinearly processing coefficients of the orthogonal or biorthog- 
onal wavelet transform, is certainly not desirable for this particular application. For 
more details on different enhancement strategies in mammography, please refer to 
[11. 22, 23, 25]. 

Another problem related to enhancement and other kinds of processing of dig- 
ital mammograms (and not unique to mammography) is that the same mammogram 
could have been easily translated and/or rotated. From obvious reasons, it is not 
desirable that the result of processing can be significantly affected by positioning 
only — one more reason for a translation and rotation-invariant transform. 

Figure 6.3 shows outputs from steerable basis filters applied to the image from 
Figure 6.2(a) at scales {1.2.4}. By linearly combining these outputs analysis along 
arbitrary direction can be performed (cf. Chapter 4). This enables, as illustrated in 
Chapter 1. a rotation-invariant processing [21]. 



64 






Figure 6.2: (a) An original mammographic image containing a spicular mass. An 
enhanced image using (b) discrete biorthogonal wavelet transform (Bior6.8) and (c) 
multiscale spline derivative-based transform. 



(>.-> 




Figure 6.3: The magnitude of filter coefficients for G%'{u} T ,u<jy) at fc = {0, \t, \tc] 
(three levels of analysis shown). 



CHAPTER 7 
CONCLUSION 



In this thesis we constructed a new transform that does not introduce artifacts 
due to translation and rotation invariance, which are inherent to traditional wavelet 
analysis. 

We extended the one-dimensional discrete dyadic wavelet transform to higher- 
order derivatives and even-order spline functions and developed an improved initial- 
ization procedure. Comparison to the originally employed initialization [28] showed 
significantly better performance of our procedure for finer scales of analysis. 

We developed several two-dimensional transforms. All of them were derived 
with the goal of eliminating artifacts due to lack of translation and rotation in- 
variance. We presented a two-dimensional discrete dyadic wavelet transform with a 
first-derivative wavelet as an extension of the one originally proposed in [28], a two- 
dimensional discrete dyadic wavelet transform with a second-derivative wavelet that 
can approximate Laplacian of a. Gaussian, and a steerable dyadic wavelet transform 
implemented in a near-perfect reconstruction filter bank which may be preferred when 
there is a need for orientation processing along equally spaced angles. We derived 
a multiscale spline derivative-based transform which uses x-y separable filters in a 
perfect reconstruction filter bank and enables fast translation and rotation-invariant 
directional analysis of images. 

We empirically showed the presence of artifacts in image fusion and enhance- 
ment applications when traditional wavelet methods were used. Here developed trans- 
forms did not exhibit similar artifacts, and, in case of fusion, yielded nice results 
without the use of postprocessing. 



66 



67 

Future work will concentrate on a more thorough exploitation of the transform 
space for the existing applications, and potential new applications (e.g., denoising). In 
particular, we will try to use the causality property of a Gaussian kernel mentioned 
in Chapter I to trace features across scales and process them according to their 
multiscale behavior. 



REFERENCES 



1] A. Aldroubi, M. Unser, and M. Eden, "Cardinal spline filters: Stability and 
convergence to the ideal sine interpolator," Signal Process., vol. 28, no. 2, pp. 
127-138, 1992. 

2] J. Babaud. A. P. Witkin, M. Baudin, and R. 0. Duda, "Uniqueness of the 
Gaussian kernel for scale-space filtering," IEEE Trans. Pattern Anal. Mach. 
Intel!., vol. 8. no. 1, pp. 26-33, 1986. 

3] C. de Boor, A Practical Guide to Splines, Springer- Verlag, New York, NY, 1978. 

4] P. J. Burt and E. H. Adelson, "The Laplacian pyramid as a compact image 
code," IEEE Trans. Commun., vol. 31, no. 4, pp. 532-540, 1983. 

5] P. J. Burt and R. J. Kolczynski, "Enhanced image capture through fusion," in 
Proc. Int. Conf. Comput. Vision, Berlin, Germany, 1993, pp. 173-182. 

6] J. Canny, "A computational approach to edge detection," IEEE Trans. Pattern 
Anal. Mach. IntelL, vol. 8. no. 6, pp. 679-698, 1986. 

7] L. J. Chipman, T. M. Orr, and L. N. Graham, "Wavelets and image fusion," 
in Proc. IEEE Int. Conf. Image Process., Washington, D.C., 1995, vol. 3, pp. 

248-251. 

8] C. K. Chui, An Introduction to Wavelets, Academic Press, San Diego, CA, 1992. 

9] C. K. Chui, Ed., Wavelets: A Tutorial in Theory and Applications, Academic 
Press. San Diego, CA, 1992. 

[10] I. Daubechies, Ten Lectures on Wavelets, SIAM, Philadelphia, PA, 1992. 

[\Y\ J. Fan and A. Laine, "Multiscale contrast enhancement and denoising in digital 
radiographs," in Wavelets in Medicine and Biology, A. Aldroubi and M. Unser, 
Eds.. CRC Press, Boca Raton, FL, 1996, pp. 163-189. 

[12] W. T. Freeman and E. H. Adelson, "The design and use of steerable filters," 
IEEE Trans. Pattern Anal. Mach. IntelL, vol. 13, no. 9, pp. 891-906, 1991. 

[13] H. Greenspan, S. Belongie, R. Goodman, P. Perona, S. Rakshit, and C. H. 
Anderson, "Overcomplete steerable pyramid filters and rotation invariance," in 
Proc. IEEE Comput. Soc. Conf. Comput. Vision Pattern Recognit., Seattle, WA, 
1994, pp. 222-228. 

[14] Y. Hel-Or and P. C. Teo, "Canonical decomposition of steerable functions," in 
Proc. IEEE Comput. Soc. Conf. Comput. Vision Pattern Recognit., San Fran- 
cisco, CA, 1996, pp. 809-816. 



69 



[15] M. Holschneider, R. Kronland-Martinet, J. Morlet, and Ph. Tchamitchian, "A 
real-time algorithm for signal analysis with the help of the wavelet transform," in 
Wavelets, Time-Frequency Methods and Phase Space, J. M. Combes, A. Gross- 
mann, and Ph. Tchamitchian, Eds., Springer- Verlag, Berlin, Germany, 1989, pp. 
286-297. 

[16] A. J. Jerri, "The Shannon sampling theorem — its various extensions and appli- 
cations: A tutorial review," Proc. IEEE, vol. 65, no. 11, pp. 1565-1596, 1977. 

[17] A. Karasaridis and E. Simoncelli, "A filter design technique for steerable pyramid 
image transforms," in Proc. IEEE Int. Conf Acoust. Speech Signal Process., 
Atlanta, GA, 1996, vol. 4, pp. 2389-2392. 

[18] I. Koren and A. Laine, "A discrete dyadic wavelet transform for multidimensional 
feature analysis," in Time-Frequency and Wavelet Transforms in Biomedical 
Engineering, M. Akay, Ed.. IEEE Press, New York, NY, 1997. 

[19] I. Koren, A. F. Laine, J. Fan, and F. J. Taylor, "Edge detection in echocardio- 
graphic image sequences by 3-d multiscale analysis," in Proc. IEEE Int. Conf. 
Mage Process., Austin, TX, 1994, vol. 1, pp. 288-292. 

[20] I. Koren. A. Laine. and F. Taylor, "Image fusion using steerable dyadic wavelet 
transform," in Proc. IEEE Int. Conf. Image Process., Washington, D.C., 1995, 
vol. 3, pp. 232-235. 





21}) I. Koren, A. Laine. F. Taylor, and M. Lewis, "Interactive wavelet processing and 
techniques applied to digital mammography," in Proc. IEEE Int. Conf. Acoust. 
Speech Signal Process., Atlanta, GA, 1996, vol, 3, pp. 1415-1418. 

] A. Laine, J. Fan, and S. Schuler, "A framework for contrast enhancement by 
dyadic wavelet analysis." in Digital Mammography, A. G. Gale, S. M. Astley, 
D. R. Dance, and A. Y. Cairns, Eds., Elsevier, Amsterdam, The Netherlands, 
1994, pp. 91-100. 

A. Laine, J. Fan, and W. Yang, "Wavelets for contrast enhancement of digital 
mammography," IEEE Eng. Med. Biol. Mag., vol. 14, no. 5, pp. 536-550, 1995. 

[24] A. Laine, I. Koren, W. Yang, and F. Taylor, "A steerable dyadic wavelet 
transform and interval wavelets for enhancement of digital mammography," in 
Wavelet Applications II. Proc. SPIE, H. H. Szu, Ed., Orlando, FL, 1995, vol. 
2/91, pp. 736-749. 

[25] A. F. Laine, S. Schuler, J. Fan, and W. Huda, "Mammographic feature en- 
hancement bv multiscale analysis," IEEE Trans. Med. Imaging, vol. 13. no. 4, 
pp. 725-740, 1994. 

[26] H. Li, B. S. Manjunath, and S. K. Mitra, "Multi-sensor image fusion using the 
wavelet transform," in Proc. IEEE Int. Conf. Image Process., Austin, TX, 1994, 
vol. 1, pp. 51-55. 

[27] S. Mallat and W. L. Hwang, "Singularity detection and processing with 
wavelets," IEEE Trans. Inf. Theory, vol. 38, no. 2, pp. 617-643, 1992. 

[28] S. Mallat and S. Zhong, "Characterization of signals from multiscale edges," 
IEEE Trans. Pattern Anal. Mach. Inteli, vol. 14, no. 7, pp. 710-732, 1992. 



71) 



[29] D. Marr and E. Hildreth, "Theory of edge detection," Proc. R. Soc. London 
Ser. B, vol. 207, no. 1167. pp. 187-217, 1980. 

[30] A. V. Oppenheim and R. W. Schafer, Discrete-Time Signal Processing, Prentice- 
Hall, Englewood Cliffs, NJ, 1989. 

[31] P. Perona. "Deformable kernels for early vision," IEEE Trans. Pattern Anal. 
Mach. Intel!., vol. 17, no. 5, pp. 488-499, 1995. 

[32] T. Ranchin, L. Wald, and M. Mangolini, "Efficient data fusion using wavelet 
transform: the case of SPOT satellite images," in Matin nmtical Imaging: 
Wavelet Applications in Signal and Image Processing. Proc. SPIE. A. F. Laine, 
Ed., San Diego, CA, 1993, vol. 2034, pp. 171-178. 

[33] 0. Rioul and P. Duhamel, "Fast algorithms for discrete and continuous wavelet 
transforms," IEEE Trans. Inf. Theory, vol. 38, no. 2, pp. 569-586, 1992. 

[34] I. J. Schoenberg, "Contributions to the problem of approximation of equidistant 
data by analytic functions," Q. Appl. Math., vol. 4. no. 1. pp. 45-99, 112-141, 
1946. ' 

[35] I. J. Schoenberg. "Cardinal interpolation and spline functions." ./. Approx. 
Theory, vol. 2, no. 2, pp. 167-206, 1969. 

[36] I. J. Schoenberg, "Notes on spline functions III: On the convergence of the 
interpolating cardinal splines as their degree tends to infinity," Isr. J. Math.. 
vol. 16. no. 1, pp. 87-93, 1973. 

[37] C. E. Shannon, "Communication in the presence of noise," Proc. IRE. vol. 37, 
no. 1, pp. 10-21, 1949. 

[38] E. P. Simoncelli and W. T. Freeman, "The steerable pyramid: A flexible archi- 
tecture for multi-scale derivative computation," in Proc. IEEE Int. Conf. Image 
Process., Washington, D.C., 1995, vol. 3, pp. 444-447. 

[39] E. P. Simoncelli, W. T. Freeman, E. H. Adelson, and D. J. Heeger, "Shiftable 
multiscale transforms," IEEE Trans. Inf. Theory, vol. 38, no. 2, pp. 587-607, 
1992. 

[40] M. J. T. Smith and T. P. Barnwell, III, "Exact reconstruction techniques for 
tree-structured subband coders," IEEE Trans. Acoust. Speech Signal Process., 
vol. 34, no. 3, pp. 434-441, 1986. 

[41] A. Toet, "Image fusion by a ratio of low-pass pyramid." Pattern Recognit. Lett., 
vol. 9. no. 4, pp. 245-253, 1989. 

[42] M. Unser and A. Aldroubi, "A general sampling theory for nonideal acquisition 
devices," IEEE Trans. Signal Process., vol. 42, no. 11, pp. 2915-2925, 1994. 

[43] M. Unser, A. Aldroubi, and M. Eden, "Fast B-spline transforms for continuous 
image representation and interpolation," IEEE Trans. Pattern Anal. Mach. 
Intel!., vol. 13. no. 3, pp. 277-285, 1991. 

[44] M. Unser, A. Aldroubi. and M. Eden, "On the asymptotic convergence of B- 
spline wavelets to Gabor functions," IEEE Trans. Inf. Theory, vol. 38, no. 2, 
pp. 864-872, 1992. 



71 



[45] M. Unser, A. Aldroubi, and M. Eden, "Polynomial spline signal approximations: 
Filter design and asymptotic equivalence with Shannon's sampling theorem," 
IEEE Trans. Inf. Theory, vol. 38, no. 1, pp. 95-103, 1992. 

[46] M. Unser, A. Aldroubi, and S. J. Schiff, "Fast implementation of the continuous 
wavelet transform with integer scales," IEEE Trans. Signal Process., vol. 42, no. 
12, pp. 3519-3523, 1994. 

[47] M. Vrhel, C. Lee, and M. Unser, "Fast computation of the continuous wavelet 
transform through oblique projections," in Proc. IEEE Int. Conf. Acoust. Speech 
Signal Process., Atlanta, GA, 1996, vol. 3, pp. 1459-1462. 

[48] A. Witkin. "Scale space filtering," in Proc. Int. Joint Conf. Artif. IntelL, 
Karlsruhe, Germany, 1983, pp. 1019-1022. 



BIOGRAPHICAL SKETCH 

Iztok Koren earned the B.S. and M.S. degrees in electrical engineering from the 
University of Ljubljana, Slovenia, in 1987 and 1991, respectively. He will receive the 
Ph.D. degree in electrical and computer engineering from the University of Florida, 
Gainesville, in 1996. 

His research interests include image processing, digital signal processing, and 
wavelets. 

He is a member of Eta Kappa Nu, Tau Beta Pi, and the IEEE. 



7-2 



I certify that I have read this study and that in my opinion it conforms to 
acceptable standards of scholarly presentation and is fully adequate, in scope and 
quality, as a dissertation for the degree of Doctor of Philosophy. 

h/4/ /f4f/fr^ 



,Fred J/Taylor/, Chairman 
Profess'or of/Electrical and Computer 
Engineering 



I certify that I have read this study and that in my opinion it conforms to 
acceptable standards of scholarly presentation and is fully adequate, in scope and 
quality, as a dissertation for the degree of Doctor of Philosophy. 



6^/iU. yj , -<^£Ui> 



Andrew F. Laine , Cochairman 
Associate Professor of Computer and 
Information Science and Engineering 



I certify that I have read this study and that in my opinion it conforms to 
acceptable standards of scholarly presentation and-^-fylly adequate, in scope and 
quality, as a dissertation for the degree of Doctor 



/ 




Jose C/ 

olr of Ele^rical and Computer 
ngmeering 



I certify that I have read this study and that in my opinion it conforms to 
acceptable standards of scholarly presentation and is fully adequate, in scope and 
quality, as a dissertation for the degree of DocJ*^ of Philosophy. 

'John M. M. Anderson 
Assistant Professor of Electrical and 
Computer Engineering 

I certify that I have read this study and that in my opinion it conforms to 
acceptable standards of scholarly presentation and is fully adequate, in scope and 
quality, as a dissertation for the degree of Doyetor of Philosophy. 






■pai 



6/lA 



?mit Sigmon 
Associate Professor of Mathematics 



This dissertation was submitted to the Graduate Faculty of the College of 
Engineering and to the Graduate School and was accepted as partial fulfillment of 
the requirements for the degree of Doctor of Philosophy. 



December 1996 



V^ Winfred M. Phillips 

Dean, College of Engineering 



Karen A. Holbrook 
Dean, Graduate School 



L 



LD 

1780 

19S£ 

.Km 



UNIVERSITY OF FLORIDA 



3 1262 08555 0696